EmoSSR: face ratings

Here we assess whether the faces used in the experiments conveyed the intended emotional expression.

# setwd(paste0(getwd(), "/ratings")) # set working subdir

# pilot1 must be treated separately because it used different faces
pilot1 <- read_csv("pilot1_ratings.csv") %>%
  mutate(
    exp = "pilot1",
    participant = as.factor(subject),
    question = as.factor(question),
    emotion = factor(substr(pic_id, 1, 1), labels = c("anger", "neutral", "happiness")),
    picID = as.factor(substr(pic_id, 2, 2))
  ) %>%
  dplyr::select(exp, participant, question, picID, emotion, ratings = "rating")

exps <- c( # variable with experiment identifier
  # participants x picID x question
  rep(sub("_ratings.csv", "", list.files(pattern = ".csv")[1]), each = (30 * 8 * 3)),
  rep(sub("_ratings.csv", "", list.files(pattern = ".csv")[2]), each = (30 * 8 * 3)),
  rep(sub("_ratings.csv", "", list.files(pattern = ".csv")[4]), each = (16 * 8 * 3))
)

ratings <- list.files(pattern = ".csv")[c(1, 2, 4)] %>% # list of csv files in current directory (excluding pilot1)
  map_df(~ read_csv(.)) %>%
  mutate(
    exp = as.factor(exps),
    participant = as.factor(subject),
    question = as.factor(question),
    picID = as.factor(pic_id),
    emotion = case_when(
      picID == "1" | picID == "5" ~ "anger",
      picID == "2" | picID == "10" ~ "neutral",
      picID == "3" | picID == "15" ~ "happiness",
      picID == "4" | picID == "20" ~ "disgust"
    )
  ) %>%
  dplyr::select(exp, participant, question, picID, emotion, ratings = "rating") %>%
  bind_rows(., pilot1) %>%
  mutate(
    variability = case_when(
      exp != "pilot1" & picID == "1" |
        exp != "pilot1" & picID == "2" |
        exp != "pilot1" & picID == "3" |
        exp != "pilot1" & picID == "4" ~ "low",
      exp != "pilot1" & picID == "5" |
        exp != "pilot1" & picID == "10" |
        exp != "pilot1" & picID == "15" |
        exp != "pilot1" & picID == "20" ~ "high",
      exp == "pilot1" ~ "NA"
    )
  ) %>%
  unite(condition, c(emotion, variability), sep = "_", remove = FALSE) %>%
  mutate(
    exp = factor(exp,
      levels = c("pilot1", "pilot2", "exp1", "exp2") # sort pilots first
    ),
    condition = recode(factor(condition),
      "anger_low" = "anger, low variability",
      "neutral_low" = "neutral, low variability",
      "happiness_low" = "happiness, low variability",
      "disgust_low" = "disgust, low variability",
      "anger_high" = "anger, high variability",
      "neutral_high" = "neutral, high variability",
      "happiness_high" = "happiness, high variability",
      "disgust_high" = "disgust, high variability",
      "anger_NA" = "anger",
      "neutral_NA" = "neutral",
      "happiness_NA" = "happiness"
    )
  ) %>%
  dplyr::select(-picID)

rm(pilot1)

Face ratings were collected at the end of all experiments. In Pilot 1, we used 9 identities and 2 emotional expressions, whereas in the other experiments we used 2 identities and 3 expressions (disgust as a filler, never regularly presented during trials).

summary.ratings <- summarySEwithin(
  data = ratings,
  measurevar = "ratings",
  betweenvars = "exp",
  withinvars = c("condition", "question"),
  idvar = "participant"
) %>%
  mutate(
    CI.95.low = ratings - ci,
    CI.95.high = ratings + ci
  ) %>%
  dplyr::select(exp, condition, question, ratings, CI.95.low, CI.95.high) %>%
  mutate(condition = factor(condition,
    levels = c(
      "neutral", "anger", "happiness",
      "neutral, low variability", "anger, low variability", "happiness, low variability", "disgust, low variability", 
      "neutral, high variability", "anger, high variability", "happiness, high variability", "disgust, high variability"
    )
  )) %>%
  arrange(exp, question, condition)

kable(summary.ratings, digits = 2)

exp	condition	question	ratings	CI.95.low	CI.95.high
pilot1	neutral	angry	4.00	3.73	4.27
pilot1	anger	angry	7.02	6.72	7.33
pilot1	happiness	angry	1.13	0.95	1.30
pilot1	neutral	happy	3.18	2.97	3.40
pilot1	anger	happy	1.44	1.20	1.68
pilot1	happiness	happy	6.86	6.62	7.11
pilot2	neutral, low variability	angry	3.02	2.17	3.88
pilot2	anger, low variability	angry	5.96	5.20	6.72
pilot2	happiness, low variability	angry	1.34	1.04	1.63
pilot2	disgust, low variability	angry	3.40	2.35	4.44
pilot2	neutral, high variability	angry	4.40	3.32	5.47
pilot2	anger, high variability	angry	8.65	8.30	9.00
pilot2	happiness, high variability	angry	1.27	1.00	1.54
pilot2	disgust, high variability	angry	4.27	2.88	5.66
pilot2	neutral, low variability	disgusted	3.09	2.07	4.10
pilot2	anger, low variability	disgusted	4.65	3.74	5.55
pilot2	happiness, low variability	disgusted	1.59	1.24	1.93
pilot2	disgust, low variability	disgusted	7.71	6.90	8.52
pilot2	neutral, high variability	disgusted	3.90	2.86	4.93
pilot2	anger, high variability	disgusted	3.59	2.52	4.65
pilot2	happiness, high variability	disgusted	1.59	1.23	1.94
pilot2	disgust, high variability	disgusted	7.71	6.47	8.95
pilot2	neutral, low variability	happy	3.52	2.66	4.39
pilot2	anger, low variability	happy	1.71	1.21	2.21
pilot2	happiness, low variability	happy	7.46	6.71	8.21
pilot2	disgust, low variability	happy	2.27	1.47	3.08
pilot2	neutral, high variability	happy	2.77	2.06	3.49
pilot2	anger, high variability	happy	1.27	0.92	1.63
pilot2	happiness, high variability	happy	7.96	7.32	8.61
pilot2	disgust, high variability	happy	1.46	0.96	1.97
exp1	neutral, low variability	angry	2.29	1.87	2.70
exp1	anger, low variability	angry	5.92	5.13	6.72
exp1	happiness, low variability	angry	1.52	1.09	1.96
exp1	disgust, low variability	angry	3.86	3.11	4.61
exp1	neutral, high variability	angry	3.59	3.05	4.12
exp1	anger, high variability	angry	8.46	7.86	9.05
exp1	happiness, high variability	angry	1.39	1.15	1.63
exp1	disgust, high variability	angry	4.79	3.82	5.76
exp1	neutral, low variability	disgusted	2.79	2.36	3.22
exp1	anger, low variability	disgusted	4.22	3.52	4.92
exp1	happiness, low variability	disgusted	2.06	1.36	2.76
exp1	disgust, low variability	disgusted	7.52	6.65	8.39
exp1	neutral, high variability	disgusted	3.39	2.78	4.00
exp1	anger, high variability	disgusted	3.62	2.86	4.38
exp1	happiness, high variability	disgusted	2.42	1.61	3.24
exp1	disgust, high variability	disgusted	7.76	6.85	8.66
exp1	neutral, low variability	happy	2.86	2.40	3.31
exp1	anger, low variability	happy	2.02	1.63	2.41
exp1	happiness, low variability	happy	8.36	7.93	8.79
exp1	disgust, low variability	happy	1.96	1.56	2.35
exp1	neutral, high variability	happy	2.16	1.77	2.55
exp1	anger, high variability	happy	1.72	1.27	2.17
exp1	happiness, high variability	happy	8.29	7.88	8.70
exp1	disgust, high variability	happy	1.59	1.33	1.85
exp2	neutral, low variability	angry	2.60	2.15	3.06
exp2	anger, low variability	angry	6.44	5.87	7.00
exp2	happiness, low variability	angry	1.27	0.91	1.63
exp2	disgust, low variability	angry	3.80	3.13	4.47
exp2	neutral, high variability	angry	4.14	3.43	4.85
exp2	anger, high variability	angry	8.50	7.97	9.03
exp2	happiness, high variability	angry	1.64	1.10	2.17
exp2	disgust, high variability	angry	3.94	3.09	4.78
exp2	neutral, low variability	disgusted	2.40	1.86	2.95
exp2	anger, low variability	disgusted	4.54	3.90	5.18
exp2	happiness, low variability	disgusted	1.57	1.26	1.87
exp2	disgust, low variability	disgusted	7.80	7.06	8.55
exp2	neutral, high variability	disgusted	3.77	2.98	4.56
exp2	anger, high variability	disgusted	3.14	2.39	3.88
exp2	happiness, high variability	disgusted	1.87	1.23	2.51
exp2	disgust, high variability	disgusted	8.07	7.39	8.74
exp2	neutral, low variability	happy	3.70	3.14	4.27
exp2	anger, low variability	happy	1.90	1.58	2.23
exp2	happiness, low variability	happy	7.87	7.22	8.52
exp2	disgust, low variability	happy	1.50	1.29	1.72
exp2	neutral, high variability	happy	3.00	2.30	3.70
exp2	anger, high variability	happy	1.50	1.05	1.96
exp2	happiness, high variability	happy	8.24	7.81	8.66
exp2	disgust, high variability	happy	1.37	1.15	1.58

ratings.plot <- 
  ggplot(
    ratings,
    aes(
      x = question,
      y = ratings,
      color = emotion,
      fill = emotion
    )
  ) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = TRUE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(
    limits = c(1, 9),
    breaks = seq(1, 9, 1)
  ) +
  geom_hline(
    yintercept = seq(1, 9, 1),
    linetype = "dotted",
    colour = "#999999",
    size = .8,
    alpha = .5
  ) +
  labs(
    x = "question",
    y = "ratings"
  ) +
  facet_wrap(~exp, scales = "free") +
  theme_EmoSSR

ratings.plot

# ### save plots to file
# # all plots in one figure
# save_plot(paste0(getwd(), "/ratings.svg"),
#           ratings.plot,
#           base_aspect_ratio = 1,
#           base_height = 11,
#           base_width = 17
#           )

For the statistical analyses, we will fit Bayesian multilevel models using brms, separately for each question. The constant effect will be emotion. Intercepts and slopes will be allowed to vary as a function of participant and picID.

Pilot 1

• 9 face identities
• 2 emotions (anger, happiness)

How ANGRY are these faces?

# subset data
pilot1.Qangry.ratings <- ratings %>%
  dplyr::select(-c(variability, emotion)) %>%
  filter(exp == "pilot1", question == "angry") %>%
  droplevels() %>%
  mutate(condition = factor(condition,
    levels = c("anger", "neutral", "happiness") # sort anger first
  ))

pilot1.Qangry.model <- brm(ratings ~ condition # constant effects
+ (condition | participant), # varying effects
data = pilot1.Qangry.ratings, # subset data
family = cumulative("probit"), # likelihood function
prior = weak.priors, # priors
sample_prior = TRUE, # draw samples from the prior distribution
inits = "0", # initial parameter values in the chains
control = list(
  adapt_delta = .99, # parameters of NUTS sampler
  max_treedepth = 15
),
chains = num.chains, # number of chains
iter = num.iter, # number of iterations
warmup = num.warmup, # number of warm-up samples
thin = num.thin, # thinning rate
algorithm = "sampling", # type of sampling algorithm
cores = num.chains, # number of processor cores to use when running chains in parallel
seed = seed.smorfia # RNG seed
)

saveRDS(pilot1.Qangry.model, file = "pilot1_Qangry_model.rds", compress = "gzip") # save results

# subset data
pilot1.Qangry.ratings <- ratings %>%
  dplyr::select(-c(variability, emotion)) %>%
  filter(exp == "pilot1", question == "angry") %>%
  droplevels() %>%
  mutate(condition = factor(condition,
    levels = c("anger", "neutral", "happiness") # sort anger first
  ))

# load model
pilot1.Qangry.model <- readRDS("pilot1_Qangry_model.rds")

# summary posterior distributions & diagnostics
summary.pilot1.Qangry.model <-
  describe_posterior(pilot1.Qangry.model,
    centrality = "mean",
    ci = 0.95,
    ci_method = "hdi",
    test = NULL,
    diagnostic = "all"
  ) %>%
  as_tibble() %>%
  separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
  dplyr::select(-coef) %>%
  mutate(Parameter = recode(factor(Parameter),
    "Intercept.1." = "Intercept[1]",
    "Intercept.2." = "Intercept[2]",
    "Intercept.3." = "Intercept[3]",
    "Intercept.4." = "Intercept[4]",
    "Intercept.5." = "Intercept[5]",
    "Intercept.6." = "Intercept[6]",
    "Intercept.7." = "Intercept[7]",
    "Intercept.8." = "Intercept[8]",
    "conditionneutral" = "anger minus neutral",
    "conditionhappiness" = "anger minus happiness"
  )) %>%
  dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)

kable(summary.pilot1.Qangry.model, digits = 2)

Parameter	Mean	HDI95_Low	HDI95_High	Estimated_Sample_Size	MonteCarlo_StdErr
Intercept[1]	-3.57	-4.02	-3.12	10162	0
Intercept[2]	-2.90	-3.32	-2.47	10259	0
Intercept[3]	-2.42	-2.83	-2.02	10291	0
Intercept[4]	-2.00	-2.39	-1.60	10680	0
Intercept[5]	-1.24	-1.61	-0.88	10587	0
Intercept[6]	-0.69	-1.04	-0.34	10987	0
Intercept[7]	-0.18	-0.51	0.16	10937	0
Intercept[8]	0.46	0.11	0.79	10941	0
anger minus neutral	-1.91	-2.37	-1.48	9933	0
anger minus happiness	-4.05	-4.68	-3.41	10251	0

# MCMC chains
pilot1.Qangry.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select(
    "Chain" = ".chain",
    "Intercept[1]" = "b_Intercept[1]",
    "Intercept[2]" = "b_Intercept[2]",
    "Intercept[3]" = "b_Intercept[3]",
    "Intercept[4]" = "b_Intercept[4]",
    "Intercept[5]" = "b_Intercept[5]",
    "Intercept[6]" = "b_Intercept[6]",
    "Intercept[7]" = "b_Intercept[7]",
    "Intercept[8]" = "b_Intercept[8]",
    "anger minus neutral" = "b_conditionneutral",
    "anger minus happiness" = "b_conditionhappiness"
  ) %>%
  as.data.frame() %>%
  mcmc_trace(facet_args = list(ncol = 2)) +
  geom_vline(
    xintercept = seq(1, 9, 1),
    linetype = "dotted",
    colour = "#999999",
    size = .8,
    alpha = .5
  ) +
  ggtitle("Chain convergence") +
  theme_EmoSSR +
  theme(legend.position = "none")

### observed vs. predicted data ###
pilot1.Qangry.ratings %>%
  mutate(predicted = posterior_predict(pilot1.Qangry.model)[1, ]) %>% # take only first row of sampled ratings (as example)
  dplyr::rename(observed = ratings) %>%
  gather(key = data_type, value = ratings, observed:predicted) %>%
  ggplot(aes(
    x = condition,
    y = ratings,
    color = condition,
    fill = condition
  )) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(
    limits = c(1, 9),
    breaks = seq(1, 9, 1)
  ) +
  geom_hline(
    yintercept = seq(1, 9, 1),
    linetype = "dotted",
    colour = "#999999",
    size = .8,
    alpha = .5
  ) +
  facet_wrap(~data_type, scales = "free") +
  xlab("") +
  ggtitle("pilot 1, \"anger\" question ") +
  theme_EmoSSR

pilot1.Qangry.posterior.test <- pilot1.Qangry.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select(
    "anger minus neutral" = "b_conditionneutral",
    "anger minus happiness" = "b_conditionhappiness"
  ) %>%
  gather(key = "condition", value = "ratings") %>%
  mutate(condition = as.factor(condition))

# summary
summary.pilot1.Qangry.posterior.test <- pilot1.Qangry.posterior.test %>%
  group_by(condition) %>%
  dplyr::summarize(
    mean = mean(ratings),
    hdi.95.low = hdi(ratings)[1],
    hdi.95.high = hdi(ratings)[2],
    evidence.ratio = ifelse(mean < 0,
      length(which(ratings < 0)) / length(which(ratings > 0)),
      length(which(ratings > 0)) / length(which(ratings < 0))
    )
  ) %>%
  ungroup()

kable(summary.pilot1.Qangry.posterior.test, digits = 2)

condition	mean	hdi.95.low	hdi.95.high	evidence.ratio
anger minus happiness	-4.05	-4.68	-3.41	Inf
anger minus neutral	-1.91	-2.37	-1.48	Inf

par(mfrow = c(1, 2))

for (i in 1:length(levels(pilot1.Qangry.posterior.test$condition))) {
  temp.data <- filter(pilot1.Qangry.posterior.test, condition == levels(condition)[i])

  plotPost(
    pull(temp.data, ratings),
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "",
    col = "#b3cde0",
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
}

Take-home message:

• happy and neutral faces are rated as less angry than angry faces

How HAPPY are these faces?

# subset data
pilot1.Qhappy.ratings <- ratings %>%
  dplyr::select(-c(variability, emotion)) %>%
  filter(exp == "pilot1", question == "happy") %>%
  droplevels() %>%
  mutate(condition = factor(condition,
    levels = c("happiness", "neutral", "anger") # sort happiness first
  ))

pilot1.Qhappy.model <- brm(ratings ~ condition
+ (condition | participant),
data = pilot1.Qhappy.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
  adapt_delta = .99,
  max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)

saveRDS(pilot1.Qhappy.model, file = "pilot1_Qhappy_model.rds", compress = "gzip")

# subset data
pilot1.Qhappy.ratings <- ratings %>%
  dplyr::select(-c(variability, emotion)) %>%
  filter(exp == "pilot1", question == "happy") %>%
  droplevels() %>%
  mutate(condition = factor(condition,
    levels = c("happiness", "neutral", "anger") # sort happiness first
  ))

# load model
pilot1.Qhappy.model <- readRDS("pilot1_Qhappy_model.rds")

# summary posterior distributions & diagnostics
summary.pilot1.Qhappy.model <-
  describe_posterior(pilot1.Qhappy.model,
    centrality = "mean",
    ci = 0.95,
    ci_method = "hdi",
    test = NULL,
    diagnostic = "all"
  ) %>%
  as_tibble() %>%
  separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
  dplyr::select(-coef) %>%
  mutate(Parameter = recode(factor(Parameter),
    "Intercept.1." = "Intercept[1]",
    "Intercept.2." = "Intercept[2]",
    "Intercept.3." = "Intercept[3]",
    "Intercept.4." = "Intercept[4]",
    "Intercept.5." = "Intercept[5]",
    "Intercept.6." = "Intercept[6]",
    "Intercept.7." = "Intercept[7]",
    "Intercept.8." = "Intercept[8]",
    "conditionneutral" = "happiness minus neutral",
    "conditionanger" = "happiness minus anger"
  )) %>%
  dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)

kable(summary.pilot1.Qhappy.model, digits = 2)

Parameter	Mean	HDI95_Low	HDI95_High	Estimated_Sample_Size	MonteCarlo_StdErr
Intercept[1]	-3.83	-4.31	-3.34	7737	0
Intercept[2]	-3.22	-3.70	-2.75	7787	0
Intercept[3]	-2.63	-3.09	-2.18	7790	0
Intercept[4]	-2.11	-2.55	-1.67	7885	0
Intercept[5]	-1.12	-1.52	-0.73	8214	0
Intercept[6]	-0.73	-1.10	-0.35	8914	0
Intercept[7]	0.01	-0.36	0.37	9142	0
Intercept[8]	0.79	0.42	1.18	9953	0
happiness minus neutral	-2.54	-2.97	-2.12	10829	0
happiness minus anger	-3.96	-4.49	-3.40	9539	0

# MCMC chains
pilot1.Qhappy.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select("Chain" = ".chain",
                "Intercept[1]" = "b_Intercept[1]",
                "Intercept[2]" = "b_Intercept[2]",
                "Intercept[3]" = "b_Intercept[3]",
                "Intercept[4]" = "b_Intercept[4]",
                "Intercept[5]" = "b_Intercept[5]",
                "Intercept[6]" = "b_Intercept[6]",
                "Intercept[7]" = "b_Intercept[7]",
                "Intercept[8]" = "b_Intercept[8]",
                "happiness minus neutral" = "b_conditionneutral",
                "happiness minus anger" = "b_conditionanger"
  ) %>%
  as.data.frame() %>%
  mcmc_trace(facet_args = list(ncol = 2)) +
  geom_vline(
    xintercept = seq(1, 9, 1),
    linetype = "dotted",
    colour = "#999999",
    size = .8,
    alpha = .5
  ) +
  ggtitle("Chain convergence") +
  theme_EmoSSR +
  theme(legend.position = "none")

### observed vs. predicted data ###
pilot1.Qhappy.ratings %>%
  mutate(predicted = posterior_predict(pilot1.Qhappy.model)[1, ]) %>% # take only first row of sampled ratings (as example)
  dplyr::rename(observed = ratings) %>%
  gather(key = data_type, value = ratings, observed:predicted) %>%
  ggplot(aes(x = condition, 
             y = ratings,
             color = condition,
             fill = condition)) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(limits = c(1, 9),
                     breaks = seq(1, 9, 1)) +
  geom_hline(yintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5) +
  facet_wrap( ~ data_type, scales = "free") +
  xlab("") +
  ggtitle("pilot 1, \"happiness\" question ") +
  theme_EmoSSR

pilot1.Qhappy.posterior.test <- pilot1.Qhappy.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select("happiness minus neutral" = "b_conditionneutral",
                "happiness minus anger" = "b_conditionanger"
  ) %>%
  gather(key = "condition", value = "ratings") %>%
  mutate(condition = as.factor(condition))

# summary
summary.pilot1.Qhappy.posterior.test <- pilot1.Qhappy.posterior.test %>%
  group_by(condition) %>%
  dplyr::summarize(
    mean = mean(ratings),
    hdi.95.low = hdi(ratings)[1],
    hdi.95.high = hdi(ratings)[2],
    evidence.ratio = ifelse(mean < 0, 
                            length(which(ratings < 0)) / length(which(ratings > 0)),
                            length(which(ratings > 0)) / length(which(ratings < 0))
    )
  ) %>%
  ungroup()

kable(summary.pilot1.Qhappy.posterior.test, digits = 2)

condition	mean	hdi.95.low	hdi.95.high	evidence.ratio
happiness minus anger	-3.96	-4.49	-3.40	Inf
happiness minus neutral	-2.54	-2.97	-2.12	Inf

par(mfrow = c(1, 2))

for(i in 1:length(levels(pilot1.Qhappy.posterior.test$condition))) {
  
  temp.data <- filter(pilot1.Qhappy.posterior.test, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

Take-home message:

• angry and neutral faces are rated as less happy than happy faces

---

---

Pilot 2

• 2 face identities
• 3 emotions (anger, happiness, disgust)

How ANGRY are these faces?

# subset data
pilot2.Qangry.ratings <- ratings %>%
  filter(exp == "pilot2", question == "angry") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
  ))

pilot2.Qangry.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = pilot2.Qangry.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
  adapt_delta = .99,
  max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)

saveRDS(pilot2.Qangry.model, file = "pilot2_Qangry_model.rds", compress = "gzip")

# subset data
pilot2.Qangry.ratings <- ratings %>%
  filter(exp == "pilot2", question == "angry") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
  ))

# load model
pilot2.Qangry.model <- readRDS("pilot2_Qangry_model.rds")

# summary posterior distributions & diagnostics
summary.pilot2.Qangry.model <-
  describe_posterior(pilot2.Qangry.model,
    centrality = "mean",
    ci = 0.95,
    ci_method = "hdi",
    test = NULL,
    diagnostic = "all"
  ) %>%
  as_tibble() %>%
  separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
  dplyr::select(-coef) %>%
  mutate(Parameter = recode(factor(Parameter),
    "Intercept.1." = "Intercept[1]",
    "Intercept.2." = "Intercept[2]",
    "Intercept.3." = "Intercept[3]",
    "Intercept.4." = "Intercept[4]",
    "Intercept.5." = "Intercept[5]",
    "Intercept.6." = "Intercept[6]",
    "Intercept.7." = "Intercept[7]",
    "Intercept.8." = "Intercept[8]",
    "emotionanger.variabilityhigh" = "anger, high variability",
    "emotionneutral.variabilityhigh" = "neutral, high variability",
    "emotionhappiness.variabilityhigh" = "happiness, high variability",
    "emotiondisgust.variabilityhigh" = "disgust, high variability",
    "emotionanger.variabilitylow" = "anger, low variability",
    "emotionneutral.variabilitylow" = "neutral, low variability",
    "emotionhappiness.variabilitylow" = "happiness, low variability",
    "emotiondisgust.variabilitylow" = "disgust, low variability"
  )) %>%
  dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)

kable(summary.pilot2.Qangry.model, digits = 2)

Parameter	Mean	HDI95_Low	HDI95_High	Estimated_Sample_Size	MonteCarlo_StdErr
Intercept[1]	-2.29	-4.79	0.04	7398	0.01
Intercept[2]	-0.65	-3.08	1.59	8242	0.01
Intercept[3]	0.19	-2.13	2.53	8451	0.01
Intercept[4]	0.83	-1.47	3.23	8392	0.01
Intercept[5]	1.45	-0.93	3.78	8145	0.01
Intercept[6]	2.79	0.27	5.22	7620	0.01
Intercept[7]	4.78	2.07	7.59	6950	0.02
Intercept[8]	7.24	4.25	10.47	7065	0.02
anger, high variability	7.44	4.34	10.38	6673	0.02
neutral, high variability	0.58	-1.80	2.84	8635	0.01
happiness, high variability	-4.55	-7.37	-1.82	8455	0.02
disgust, high variability	0.16	-2.41	2.76	8806	0.01
anger, low variability	2.36	0.05	4.83	7834	0.01
neutral, low variability	-0.99	-3.32	1.24	7900	0.01
happiness, low variability	-4.68	-7.69	-1.76	8966	0.02
disgust, low variability	-0.34	-2.79	2.04	9104	0.01

# MCMC chains
pilot2.Qangry.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select("Chain" = ".chain",
                "Intercept[1]" = "b_Intercept[1]",
                "Intercept[2]" = "b_Intercept[2]",
                "Intercept[3]" = "b_Intercept[3]",
                "Intercept[4]" = "b_Intercept[4]",
                "Intercept[5]" = "b_Intercept[5]",
                "Intercept[6]" = "b_Intercept[6]",
                "Intercept[7]" = "b_Intercept[7]",
                "Intercept[8]" = "b_Intercept[8]",
                "anger, high variability" = "b_emotionanger:variabilityhigh",
                "neutral, high variability" = "b_emotionneutral:variabilityhigh",
                "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
                "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
                "anger, low variability" = "b_emotionanger:variabilitylow",
                "neutral, low variability" = "b_emotionneutral:variabilitylow",
                "happiness, low variability" = "b_emotionhappiness:variabilitylow",
                "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  as.data.frame() %>%
  mcmc_trace(facet_args = list(ncol = 2)) +
  geom_vline(xintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5
  ) +
  ggtitle("Chain convergence") +
  theme_EmoSSR +
  theme(legend.position = "none")

### observed vs. predicted data ###
pilot2.Qangry.ratings %>%
  dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
  mutate(predicted = posterior_predict(pilot2.Qangry.model)[1, ]) %>% # take only first row of sampled ratings (as example)
  dplyr::rename(observed = ratings) %>%
  gather(key = data_type, value = ratings, observed:predicted) %>%
  ggplot(aes(x = condition, 
             y = ratings,
             color = condition,
             fill = condition)) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(limits = c(1, 9),
                     breaks = seq(1, 9, 1)) +
  geom_hline(yintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5) +
  facet_wrap( ~ data_type, scales = "free") +
  xlab("") +
  ggtitle("pilot 2, \"anger\" question ") +
  theme_EmoSSR

pilot2.Qangry.posterior.test <- pilot2.Qangry.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select(
    "anger, high variability" = "b_emotionanger:variabilityhigh",
    "neutral, high variability" = "b_emotionneutral:variabilityhigh",
    "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
    "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
    "anger, low variability" = "b_emotionanger:variabilitylow",
    "neutral, low variability" = "b_emotionneutral:variabilitylow",
    "happiness, low variability" = "b_emotionhappiness:variabilitylow",
    "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  mutate(
    ### emotion (anger minus others)
    ## high variability
    # anger minus neutral
    diff.high.angerVSneutral = `anger, high variability` - `neutral, high variability`,
    # anger minus happiness
    diff.high.angerVShappiness = `anger, high variability` - `happiness, high variability`,
    # anger minus disgust
    diff.high.angerVSdisgust = `anger, high variability` - `disgust, high variability`,
    ## low variability
    # anger minus neutral
    diff.low.angerVSneutral = `anger, low variability` - `neutral, low variability`,
    # anger minus happiness
    diff.low.angerVShappiness = `anger, low variability` - `happiness, low variability`,
    # anger minus disgust
    diff.low.angerVSdisgust = `anger, low variability` - `disgust, low variability`,
    ### variability (high minus low)
    # anger
    diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
    # neutral
    diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
    # happiness
    diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
    # disgust
    diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
  ) %>%
  gather(key = "condition", value = "ratings") %>%
  filter(condition == c(
    "diff.high.angerVSneutral",
    "diff.high.angerVShappiness",
    "diff.high.angerVSdisgust",
    "diff.low.angerVSneutral",
    "diff.low.angerVShappiness",
    "diff.low.angerVSdisgust",
    "diff.anger.highVSlow",
    "diff.neutral.highVSlow",
    "diff.happiness.highVSlow",
    "diff.disgust.highVSlow"
  )) %>%
  mutate(condition = factor(condition,
    levels = c(
      "diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
      "diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust",
      "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
    )
  ))

# summary
summary.pilot2.Qangry.posterior.test <- pilot2.Qangry.posterior.test %>%
  group_by(condition) %>%
  dplyr::summarize(
    mean = mean(ratings),
    hdi.95.low = hdi(ratings)[1],
    hdi.95.high = hdi(ratings)[2],
    evidence.ratio = ifelse(mean < 0, 
                            length(which(ratings < 0)) / length(which(ratings > 0)),
                            length(which(ratings > 0)) / length(which(ratings < 0))
    )
  ) %>%
  ungroup()

kable(summary.pilot2.Qangry.posterior.test, digits = 2)

condition	mean	hdi.95.low	hdi.95.high	evidence.ratio
diff.high.angerVSneutral	6.86	4.55	9.27	Inf
diff.high.angerVShappiness	12.01	8.53	15.63	Inf
diff.high.angerVSdisgust	7.30	4.60	10.14	Inf
diff.low.angerVSneutral	3.36	1.69	4.98	Inf
diff.low.angerVShappiness	7.07	4.47	9.90	Inf
diff.low.angerVSdisgust	2.71	0.96	4.66	399.00
diff.anger.highVSlow	5.10	3.28	7.00	Inf
diff.neutral.highVSlow	1.56	0.24	2.80	113.29
diff.happiness.highVSlow	0.13	-2.55	2.86	1.09
diff.disgust.highVSlow	0.52	-1.35	2.27	2.77

### Anger vs. others ###

pilot2.Qangry.posterior.test.emoref <- filter(
  pilot2.Qangry.posterior.test,
  condition == c(
    "diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
    "diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust"
  )
) %>%
  droplevels()

par(mfrow = c(3, 2))

for(i in 1:length(levels(pilot2.Qangry.posterior.test.emoref$condition))) {
  
  temp.data <- filter(pilot2.Qangry.posterior.test.emoref, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

### High vs. low variability ###

pilot2.Qangry.posterior.test.var <- filter(
  pilot2.Qangry.posterior.test,
  condition == c(
    "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
  )
) %>%
  droplevels()

par(mfrow = c(2, 2))

for(i in 1:length(levels(pilot2.Qangry.posterior.test.var$condition))) {
  
  temp.data <- filter(pilot2.Qangry.posterior.test.var, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

Take-home message:

• all faces are rated as less angry than angry faces, in both within-emotion variabilities
• angry and neutral faces are rated as angrier in high compared to low within-emotion variability
• happy and disgusted faces are rated as similarly angry in high and low within-emotion variability

---

---

How HAPPY are these faces?

# subset data
pilot2.Qhappy.ratings <- ratings %>%
  filter(exp == "pilot2", question == "happy") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
  ))

pilot2.Qhappy.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = pilot2.Qhappy.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
  adapt_delta = .99,
  max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)

saveRDS(pilot2.Qhappy.model, file = "pilot2_Qhappy_model.rds", compress = "gzip")

# subset data
pilot2.Qhappy.ratings <- ratings %>%
  filter(exp == "pilot2", question == "happy") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
  ))

# load model
pilot2.Qhappy.model <- readRDS("pilot2_Qhappy_model.rds")

# summary posterior distributions & diagnostics
summary.pilot2.Qhappy.model <-
  describe_posterior(pilot2.Qhappy.model,
    centrality = "mean",
    ci = 0.95,
    ci_method = "hdi",
    test = NULL,
    diagnostic = "all"
  ) %>%
  as_tibble() %>%
  separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
  dplyr::select(-coef) %>%
  mutate(Parameter = recode(factor(Parameter),
    "Intercept.1." = "Intercept[1]",
    "Intercept.2." = "Intercept[2]",
    "Intercept.3." = "Intercept[3]",
    "Intercept.4." = "Intercept[4]",
    "Intercept.5." = "Intercept[5]",
    "Intercept.6." = "Intercept[6]",
    "Intercept.7." = "Intercept[7]",
    "Intercept.8." = "Intercept[8]",
    "emotionanger.variabilityhigh" = "anger, high variability",
    "emotionneutral.variabilityhigh" = "neutral, high variability",
    "emotionhappiness.variabilityhigh" = "happiness, high variability",
    "emotiondisgust.variabilityhigh" = "disgust, high variability",
    "emotionanger.variabilitylow" = "anger, low variability",
    "emotionneutral.variabilitylow" = "neutral, low variability",
    "emotionhappiness.variabilitylow" = "happiness, low variability",
    "emotiondisgust.variabilitylow" = "disgust, low variability"
  )) %>%
  dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)

kable(summary.pilot2.Qhappy.model, digits = 2)

Parameter	Mean	HDI95_Low	HDI95_High	Estimated_Sample_Size	MonteCarlo_StdErr
Intercept[1]	-1.89	-4.15	0.44	8717	0.01
Intercept[2]	-0.49	-2.79	1.74	9417	0.01
Intercept[3]	0.69	-1.60	3.03	9774	0.01
Intercept[4]	1.19	-1.16	3.50	9573	0.01
Intercept[5]	2.39	-0.06	4.86	9343	0.01
Intercept[6]	3.57	0.95	6.08	8820	0.01
Intercept[7]	5.45	2.53	8.31	7592	0.02
Intercept[8]	7.34	3.90	10.94	6094	0.02
happiness, high variability	6.00	3.17	9.15	7681	0.02
neutral, high variability	-0.61	-2.92	1.63	9815	0.01
anger, high variability	-3.30	-5.67	-0.88	8862	0.01
disgust, high variability	-3.66	-6.67	-0.89	9665	0.01
happiness, low variability	5.11	2.27	7.89	8506	0.02
neutral, low variability	0.20	-2.16	2.43	9573	0.01
anger, low variability	-2.16	-4.53	0.15	9134	0.01
disgust, low variability	-1.64	-4.12	0.86	9796	0.01

# MCMC chains
pilot2.Qhappy.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select("Chain" = ".chain",
                "Intercept[1]" = "b_Intercept[1]",
                "Intercept[2]" = "b_Intercept[2]",
                "Intercept[3]" = "b_Intercept[3]",
                "Intercept[4]" = "b_Intercept[4]",
                "Intercept[5]" = "b_Intercept[5]",
                "Intercept[6]" = "b_Intercept[6]",
                "Intercept[7]" = "b_Intercept[7]",
                "Intercept[8]" = "b_Intercept[8]",
                "anger, high variability" = "b_emotionanger:variabilityhigh",
                "neutral, high variability" = "b_emotionneutral:variabilityhigh",
                "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
                "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
                "anger, low variability" = "b_emotionanger:variabilitylow",
                "neutral, low variability" = "b_emotionneutral:variabilitylow",
                "happiness, low variability" = "b_emotionhappiness:variabilitylow",
                "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  as.data.frame() %>%
  mcmc_trace(facet_args = list(ncol = 2)) +
  geom_vline(xintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5
  ) +
  ggtitle("Chain convergence") +
  theme_EmoSSR +
  theme(legend.position = "none")

### observed vs. predicted data ###
pilot2.Qhappy.ratings %>%
  dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
  mutate(predicted = posterior_predict(pilot2.Qhappy.model)[1, ]) %>% # take only first row of sampled ratings (as example)
  dplyr::rename(observed = ratings) %>%
  gather(key = data_type, value = ratings, observed:predicted) %>%
  ggplot(aes(x = condition, 
             y = ratings,
             color = condition,
             fill = condition)) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(limits = c(1, 9),
                     breaks = seq(1, 9, 1)) +
  geom_hline(yintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5) +
  facet_wrap( ~ data_type, scales = "free") +
  xlab("") +
  ggtitle("pilot 2, \"happiness\" question ") +
  theme_EmoSSR

pilot2.Qhappy.posterior.test <- pilot2.Qhappy.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select(
    "anger, high variability" = "b_emotionanger:variabilityhigh",
    "neutral, high variability" = "b_emotionneutral:variabilityhigh",
    "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
    "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
    "anger, low variability" = "b_emotionanger:variabilitylow",
    "neutral, low variability" = "b_emotionneutral:variabilitylow",
    "happiness, low variability" = "b_emotionhappiness:variabilitylow",
    "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  mutate(
    ### emotion (happiness minus others)
    ## high variability
    # happiness minus neutral
    diff.high.happinessVSneutral = `happiness, high variability` - `neutral, high variability`,
    # happiness minus anger
    diff.high.happinessVSanger = `happiness, high variability` - `anger, high variability`,
    # happiness minus disgust
    diff.high.happinessVSdisgust = `happiness, high variability` - `disgust, high variability`,
    ## low variability
    # happiness minus neutral
    diff.low.happinessVSneutral = `happiness, low variability` - `neutral, low variability`,
    # happiness minus anger
    diff.low.happinessVSanger = `happiness, low variability` - `anger, low variability`,
    # happiness minus disgust
    diff.low.happinessVSdisgust = `happiness, low variability` - `disgust, low variability`,
    ### variability (high minus low)
    # anger
    diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
    # neutral
    diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
    # happiness
    diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
    # disgust
    diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
  ) %>%
  gather(key = "condition", value = "ratings") %>%
  filter(condition == c(
    "diff.high.happinessVSneutral",
    "diff.high.happinessVSanger",
    "diff.high.happinessVSdisgust",
    "diff.low.happinessVSneutral",
    "diff.low.happinessVSanger",
    "diff.low.happinessVSdisgust",
    "diff.anger.highVSlow",
    "diff.neutral.highVSlow",
    "diff.happiness.highVSlow",
    "diff.disgust.highVSlow"
  )) %>%
  mutate(condition = factor(condition,
    levels = c(
      "diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
      "diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust",
      "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
    )
  ))

# summary
summary.pilot2.Qhappy.posterior.test <- pilot2.Qhappy.posterior.test %>%
  group_by(condition) %>%
  dplyr::summarize(
    mean = mean(ratings),
    hdi.95.low = hdi(ratings)[1],
    hdi.95.high = hdi(ratings)[2],
    evidence.ratio = ifelse(mean < 0, 
                            length(which(ratings < 0)) / length(which(ratings > 0)),
                            length(which(ratings > 0)) / length(which(ratings < 0))
    )
  ) %>%
  ungroup()

kable(summary.pilot2.Qhappy.posterior.test, digits = 2)

condition	mean	hdi.95.low	hdi.95.high	evidence.ratio
diff.high.happinessVSneutral	6.59	4.09	9.26	Inf
diff.high.happinessVSanger	9.27	6.50	12.48	Inf
diff.high.happinessVSdisgust	9.65	6.46	13.26	Inf
diff.low.happinessVSneutral	4.95	3.08	7.23	Inf
diff.low.happinessVSanger	7.31	4.89	9.96	Inf
diff.low.happinessVSdisgust	6.75	4.30	9.53	Inf
diff.anger.highVSlow	-1.13	-2.55	0.12	19.00
diff.neutral.highVSlow	-0.78	-1.91	0.37	12.22
diff.happiness.highVSlow	0.89	-0.80	2.72	6.77
diff.disgust.highVSlow	-2.01	-5.05	0.14	21.86

### Happiness vs. others ###

pilot2.Qhappy.posterior.test.emoref <- filter(
  pilot2.Qhappy.posterior.test,
  condition == c(
    "diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
    "diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust"
  )
) %>%
  droplevels()

par(mfrow = c(3, 2))

for(i in 1:length(levels(pilot2.Qhappy.posterior.test.emoref$condition))) {
  
  temp.data <- filter(pilot2.Qhappy.posterior.test.emoref, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

### High vs. low variability ###

pilot2.Qhappy.posterior.test.var <- filter(
  pilot2.Qhappy.posterior.test,
  condition == c(
    "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
  )
) %>%
  droplevels()

par(mfrow = c(2, 2))

for(i in 1:length(levels(pilot2.Qhappy.posterior.test.var$condition))) {
  
  temp.data <- filter(pilot2.Qhappy.posterior.test.var, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

Take-home message:

• all faces are rated as less happy than happy faces, in both within-emotion variabilities
• angry, neutral, and disgusted faces are rated as happier in low compared to high within-emotion variability
• happy faces are rated as similarly happy in high and low within-emotion variability

---

---

How DISGUSTED are these faces?

# subset data
pilot2.Qdisgust.ratings <- ratings %>%
  filter(exp == "pilot2", question == "disgusted") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
  ))

pilot2.Qdisgust.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = pilot2.Qdisgust.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
  adapt_delta = .99,
  max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)

saveRDS(pilot2.Qdisgust.model, file = "pilot2_Qdisgust_model.rds", compress = "gzip")

# subset data
pilot2.Qdisgust.ratings <- ratings %>%
  filter(exp == "pilot2", question == "disgusted") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
  ))

# load model
pilot2.Qdisgust.model <- readRDS("pilot2_Qdisgust_model.rds")

# summary posterior distributions & diagnostics
summary.pilot2.Qdisgust.model <-
  describe_posterior(pilot2.Qdisgust.model,
    centrality = "mean",
    ci = 0.95,
    ci_method = "hdi",
    test = NULL,
    diagnostic = "all"
  ) %>%
  as_tibble() %>%
  separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
  dplyr::select(-coef) %>%
  mutate(Parameter = recode(factor(Parameter),
    "Intercept.1." = "Intercept[1]",
    "Intercept.2." = "Intercept[2]",
    "Intercept.3." = "Intercept[3]",
    "Intercept.4." = "Intercept[4]",
    "Intercept.5." = "Intercept[5]",
    "Intercept.6." = "Intercept[6]",
    "Intercept.7." = "Intercept[7]",
    "Intercept.8." = "Intercept[8]",
    "emotionanger.variabilityhigh" = "anger, high variability",
    "emotionneutral.variabilityhigh" = "neutral, high variability",
    "emotionhappiness.variabilityhigh" = "happiness, high variability",
    "emotiondisgust.variabilityhigh" = "disgust, high variability",
    "emotionanger.variabilitylow" = "anger, low variability",
    "emotionneutral.variabilitylow" = "neutral, low variability",
    "emotionhappiness.variabilitylow" = "happiness, low variability",
    "emotiondisgust.variabilitylow" = "disgust, low variability"
  )) %>%
  dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)

kable(summary.pilot2.Qdisgust.model, digits = 2)

Parameter	Mean	HDI95_Low	HDI95_High	Estimated_Sample_Size	MonteCarlo_StdErr
Intercept[1]	-2.08	-4.33	0.15	6373	0.01
Intercept[2]	-0.95	-3.22	1.18	6575	0.01
Intercept[3]	-0.36	-2.57	1.83	6661	0.01
Intercept[4]	0.14	-2.11	2.29	6728	0.01
Intercept[5]	0.85	-1.28	3.14	6702	0.01
Intercept[6]	1.41	-0.85	3.58	6743	0.01
Intercept[7]	2.98	0.56	5.28	6567	0.01
Intercept[8]	4.92	2.05	7.64	6530	0.02
disgust, high variability	4.28	1.12	7.27	8458	0.02
neutral, high variability	-0.37	-2.67	1.77	7071	0.01
anger, high variability	-0.61	-2.81	1.55	6771	0.01
happiness, high variability	-3.06	-5.55	-0.76	6634	0.02
disgust, low variability	3.47	0.92	5.99	6854	0.02
neutral, low variability	-1.29	-3.55	1.03	6857	0.01
anger, low variability	0.30	-1.92	2.45	6670	0.01
happiness, low variability	-2.74	-5.18	-0.44	6657	0.01

# MCMC chains
pilot2.Qdisgust.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select("Chain" = ".chain",
                "Intercept[1]" = "b_Intercept[1]",
                "Intercept[2]" = "b_Intercept[2]",
                "Intercept[3]" = "b_Intercept[3]",
                "Intercept[4]" = "b_Intercept[4]",
                "Intercept[5]" = "b_Intercept[5]",
                "Intercept[6]" = "b_Intercept[6]",
                "Intercept[7]" = "b_Intercept[7]",
                "Intercept[8]" = "b_Intercept[8]",
                "anger, high variability" = "b_emotionanger:variabilityhigh",
                "neutral, high variability" = "b_emotionneutral:variabilityhigh",
                "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
                "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
                "anger, low variability" = "b_emotionanger:variabilitylow",
                "neutral, low variability" = "b_emotionneutral:variabilitylow",
                "happiness, low variability" = "b_emotionhappiness:variabilitylow",
                "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  as.data.frame() %>%
  mcmc_trace(facet_args = list(ncol = 2)) +
  geom_vline(xintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5
  ) +
  ggtitle("Chain convergence") +
  theme_EmoSSR +
  theme(legend.position = "none")

### observed vs. predicted data ###
pilot2.Qdisgust.ratings %>%
  dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
  mutate(predicted = posterior_predict(pilot2.Qdisgust.model)[1, ]) %>% # take only first row of sampled ratings (as example)
  dplyr::rename(observed = ratings) %>%
  gather(key = data_type, value = ratings, observed:predicted) %>%
  ggplot(aes(x = condition, 
             y = ratings,
             color = condition,
             fill = condition)) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(limits = c(1, 9),
                     breaks = seq(1, 9, 1)) +
  geom_hline(yintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5) +
  facet_wrap( ~ data_type, scales = "free") +
  xlab("") +
  ggtitle("pilot 2, \"disgusted\" question ") +
  theme_EmoSSR

pilot2.Qdisgust.posterior.test <- pilot2.Qdisgust.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select(
    "anger, high variability" = "b_emotionanger:variabilityhigh",
    "neutral, high variability" = "b_emotionneutral:variabilityhigh",
    "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
    "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
    "anger, low variability" = "b_emotionanger:variabilitylow",
    "neutral, low variability" = "b_emotionneutral:variabilitylow",
    "happiness, low variability" = "b_emotionhappiness:variabilitylow",
    "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  mutate(
    ### emotion (disgust minus others)
    ## high variability
    # disgust minus neutral
    diff.high.disgustVSneutral = `disgust, high variability` - `neutral, high variability`,
    # disgust minus anger
    diff.high.disgustVSanger = `disgust, high variability` - `anger, high variability`,
    # disgust minus happiness
    diff.high.disgustVShappiness = `disgust, high variability` - `happiness, high variability`,
    ## low variability
    # disgust minus neutral
    diff.low.disgustVSneutral = `disgust, low variability` - `neutral, low variability`,
    # disgust minus anger
    diff.low.disgustVSanger = `disgust, low variability` - `anger, low variability`,
    # disgust minus happiness
    diff.low.disgustVShappiness = `disgust, low variability` - `happiness, low variability`,
    ### variability (high minus low)
    # anger
    diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
    # neutral
    diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
    # happiness
    diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
    # disgust
    diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
  ) %>%
  gather(key = "condition", value = "ratings") %>%
  filter(condition == c(
    "diff.high.disgustVSneutral",
    "diff.high.disgustVSanger",
    "diff.high.disgustVShappiness",
    "diff.low.disgustVSneutral",
    "diff.low.disgustVSanger",
    "diff.low.disgustVShappiness",
    "diff.anger.highVSlow",
    "diff.neutral.highVSlow",
    "diff.happiness.highVSlow",
    "diff.disgust.highVSlow"
  )) %>%
  mutate(condition = factor(condition,
    levels = c(
      "diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
      "diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness",
      "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
    )
  ))

# summary
summary.pilot2.Qdisgust.posterior.test <- pilot2.Qdisgust.posterior.test %>%
  group_by(condition) %>%
  dplyr::summarize(
    mean = mean(ratings),
    hdi.95.low = hdi(ratings)[1],
    hdi.95.high = hdi(ratings)[2],
    evidence.ratio = ifelse(mean < 0, 
                            length(which(ratings < 0)) / length(which(ratings > 0)),
                            length(which(ratings > 0)) / length(which(ratings < 0))
    )
  ) %>%
  ungroup()

kable(summary.pilot2.Qdisgust.posterior.test, digits = 2)

condition	mean	hdi.95.low	hdi.95.high	evidence.ratio
diff.high.disgustVSneutral	4.66	2.20	7.56	799.00
diff.high.disgustVSanger	4.89	2.37	7.70	Inf
diff.high.disgustVShappiness	7.40	4.53	10.31	Inf
diff.low.disgustVSneutral	4.76	2.71	6.85	Inf
diff.low.disgustVSanger	3.19	1.65	4.81	1599.00
diff.low.disgustVShappiness	6.26	4.19	8.40	Inf
diff.anger.highVSlow	-0.90	-1.77	-0.09	44.71
diff.neutral.highVSlow	0.91	-0.30	2.16	17.60
diff.happiness.highVSlow	-0.33	-1.92	1.26	1.93
diff.disgust.highVSlow	0.81	-1.67	3.48	2.90

### Disgust vs. others ###

pilot2.Qdisgust.posterior.test.emoref <- filter(
  pilot2.Qdisgust.posterior.test,
  condition == c(
    "diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
    "diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness"
  )
) %>%
  droplevels()

par(mfrow = c(3, 2))

for(i in 1:length(levels(pilot2.Qdisgust.posterior.test.emoref$condition))) {
  
  temp.data <- filter(pilot2.Qdisgust.posterior.test.emoref, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

### High vs. low variability ###

pilot2.Qdisgust.posterior.test.var <- filter(
  pilot2.Qdisgust.posterior.test,
  condition == c(
    "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
  )
) %>%
  droplevels()

par(mfrow = c(2, 2))

for(i in 1:length(levels(pilot2.Qdisgust.posterior.test.var$condition))) {
  
  temp.data <- filter(pilot2.Qdisgust.posterior.test.var, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

Take-home message:

• all faces are rated as less disgusted than disgusted faces, in both within-emotion variabilities
• angry faces are rated as more disgusted in low compared to high within-emotion variability
• neutral faces are rated as more disgusted in high compared to low within-emotion variability
• disgusted faces are rated as similarly disgusted in high and low within-emotion variability

---

---

Experiment 1

• 2 face identities
• 3 emotions (anger, happiness, disgust)

How ANGRY are these faces?

# subset data
exp1.Qangry.ratings <- ratings %>%
  filter(exp == "exp1", question == "angry") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
  ))

exp1.Qangry.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp1.Qangry.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
  adapt_delta = .99,
  max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)

saveRDS(exp1.Qangry.model, file = "exp1_Qangry_model.rds", compress = "gzip")

# subset data
exp1.Qangry.ratings <- ratings %>%
  filter(exp == "exp1", question == "angry") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
  ))

# load model
exp1.Qangry.model <- readRDS("exp1_Qangry_model.rds")

# summary posterior distributions & diagnostics
summary.exp1.Qangry.model <-
  describe_posterior(exp1.Qangry.model,
    centrality = "mean",
    ci = 0.95,
    ci_method = "hdi",
    test = NULL,
    diagnostic = "all"
  ) %>%
  as_tibble() %>%
  separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
  dplyr::select(-coef) %>%
  mutate(Parameter = recode(factor(Parameter),
    "Intercept.1." = "Intercept[1]",
    "Intercept.2." = "Intercept[2]",
    "Intercept.3." = "Intercept[3]",
    "Intercept.4." = "Intercept[4]",
    "Intercept.5." = "Intercept[5]",
    "Intercept.6." = "Intercept[6]",
    "Intercept.7." = "Intercept[7]",
    "Intercept.8." = "Intercept[8]",
    "emotionanger.variabilityhigh" = "anger, high variability",
    "emotionneutral.variabilityhigh" = "neutral, high variability",
    "emotionhappiness.variabilityhigh" = "happiness, high variability",
    "emotiondisgust.variabilityhigh" = "disgust, high variability",
    "emotionanger.variabilitylow" = "anger, low variability",
    "emotionneutral.variabilitylow" = "neutral, low variability",
    "emotionhappiness.variabilitylow" = "happiness, low variability",
    "emotiondisgust.variabilitylow" = "disgust, low variability"
  )) %>%
  dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)

kable(summary.exp1.Qangry.model, digits = 2)

Parameter	Mean	HDI95_Low	HDI95_High	Estimated_Sample_Size	MonteCarlo_StdErr
Intercept[1]	-0.39	-2.59	1.79	9138	0.01
Intercept[2]	0.10	-2.09	2.26	9118	0.01
Intercept[3]	0.75	-1.41	2.97	9129	0.01
Intercept[4]	1.29	-0.90	3.50	8975	0.01
Intercept[5]	1.85	-0.42	4.01	8797	0.01
Intercept[6]	2.36	0.11	4.58	8603	0.01
Intercept[7]	3.19	0.95	5.51	8169	0.01
Intercept[8]	4.20	1.86	6.58	7656	0.01
anger, high variability	4.27	1.96	6.61	7643	0.01
neutral, high variability	0.47	-1.64	2.72	9162	0.01
happiness, high variability	-3.06	-5.98	-0.32	10597	0.01
disgust, high variability	1.00	-1.21	3.33	9176	0.01
anger, low variability	1.99	-0.34	4.14	8742	0.01
neutral, low variability	-0.71	-3.03	1.40	9099	0.01
happiness, low variability	-4.43	-7.70	-1.19	11415	0.02
disgust, low variability	0.48	-1.69	2.70	9327	0.01

# MCMC chains
exp1.Qangry.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select("Chain" = ".chain",
                "Intercept[1]" = "b_Intercept[1]",
                "Intercept[2]" = "b_Intercept[2]",
                "Intercept[3]" = "b_Intercept[3]",
                "Intercept[4]" = "b_Intercept[4]",
                "Intercept[5]" = "b_Intercept[5]",
                "Intercept[6]" = "b_Intercept[6]",
                "Intercept[7]" = "b_Intercept[7]",
                "Intercept[8]" = "b_Intercept[8]",
                "anger, high variability" = "b_emotionanger:variabilityhigh",
                "neutral, high variability" = "b_emotionneutral:variabilityhigh",
                "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
                "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
                "anger, low variability" = "b_emotionanger:variabilitylow",
                "neutral, low variability" = "b_emotionneutral:variabilitylow",
                "happiness, low variability" = "b_emotionhappiness:variabilitylow",
                "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  as.data.frame() %>%
  mcmc_trace(facet_args = list(ncol = 2)) +
  geom_vline(xintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5
  ) +
  ggtitle("Chain convergence") +
  theme_EmoSSR +
  theme(legend.position = "none")

### observed vs. predicted data ###
exp1.Qangry.ratings %>%
  dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
  mutate(predicted = posterior_predict(exp1.Qangry.model)[1, ]) %>% # take only first row of sampled ratings (as example)
  dplyr::rename(observed = ratings) %>%
  gather(key = data_type, value = ratings, observed:predicted) %>%
  ggplot(aes(x = condition, 
             y = ratings,
             color = condition,
             fill = condition)) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(limits = c(1, 9),
                     breaks = seq(1, 9, 1)) +
  geom_hline(yintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5) +
  facet_wrap( ~ data_type, scales = "free") +
  xlab("") +
  ggtitle("pilot 2, \"anger\" question ") +
  theme_EmoSSR

exp1.Qangry.posterior.test <- exp1.Qangry.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select(
    "anger, high variability" = "b_emotionanger:variabilityhigh",
    "neutral, high variability" = "b_emotionneutral:variabilityhigh",
    "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
    "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
    "anger, low variability" = "b_emotionanger:variabilitylow",
    "neutral, low variability" = "b_emotionneutral:variabilitylow",
    "happiness, low variability" = "b_emotionhappiness:variabilitylow",
    "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  mutate(
    ### emotion (anger minus others)
    ## high variability
    # anger minus neutral
    diff.high.angerVSneutral = `anger, high variability` - `neutral, high variability`,
    # anger minus happiness
    diff.high.angerVShappiness = `anger, high variability` - `happiness, high variability`,
    # anger minus disgust
    diff.high.angerVSdisgust = `anger, high variability` - `disgust, high variability`,
    ## low variability
    # anger minus neutral
    diff.low.angerVSneutral = `anger, low variability` - `neutral, low variability`,
    # anger minus happiness
    diff.low.angerVShappiness = `anger, low variability` - `happiness, low variability`,
    # anger minus disgust
    diff.low.angerVSdisgust = `anger, low variability` - `disgust, low variability`,
    ### variability (high minus low)
    # anger
    diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
    # neutral
    diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
    # happiness
    diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
    # disgust
    diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
  ) %>%
  gather(key = "condition", value = "ratings") %>%
  filter(condition == c(
    "diff.high.angerVSneutral",
    "diff.high.angerVShappiness",
    "diff.high.angerVSdisgust",
    "diff.low.angerVSneutral",
    "diff.low.angerVShappiness",
    "diff.low.angerVSdisgust",
    "diff.anger.highVSlow",
    "diff.neutral.highVSlow",
    "diff.happiness.highVSlow",
    "diff.disgust.highVSlow"
  )) %>%
  mutate(condition = factor(condition,
    levels = c(
      "diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
      "diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust",
      "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
    )
  ))

# summary
summary.exp1.Qangry.posterior.test <- exp1.Qangry.posterior.test %>%
  group_by(condition) %>%
  dplyr::summarize(
    mean = mean(ratings),
    hdi.95.low = hdi(ratings)[1],
    hdi.95.high = hdi(ratings)[2],
    evidence.ratio = ifelse(mean < 0, 
                            length(which(ratings < 0)) / length(which(ratings > 0)),
                            length(which(ratings > 0)) / length(which(ratings < 0))
    )
  ) %>%
  ungroup()

kable(summary.exp1.Qangry.posterior.test, digits = 2)

condition	mean	hdi.95.low	hdi.95.high	evidence.ratio
diff.high.angerVSneutral	3.79	2.86	4.79	Inf
diff.high.angerVShappiness	7.31	5.21	9.66	Inf
diff.high.angerVSdisgust	3.28	2.13	4.40	Inf
diff.low.angerVSneutral	2.73	1.73	3.69	Inf
diff.low.angerVShappiness	6.46	3.72	9.70	Inf
diff.low.angerVSdisgust	1.52	0.74	2.49	Inf
diff.anger.highVSlow	2.26	1.51	3.12	Inf
diff.neutral.highVSlow	1.18	0.39	1.92	Inf
diff.happiness.highVSlow	1.51	-1.70	5.30	4.63
diff.disgust.highVSlow	0.53	-0.39	1.46	6.11

### Anger vs. others ###

exp1.Qangry.posterior.test.emoref <- filter(
  exp1.Qangry.posterior.test,
  condition == c(
    "diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
    "diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust"
  )
) %>%
  droplevels()

par(mfrow = c(3, 2))

for(i in 1:length(levels(exp1.Qangry.posterior.test.emoref$condition))) {
  
  temp.data <- filter(exp1.Qangry.posterior.test.emoref, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

### High vs. low variability ###

exp1.Qangry.posterior.test.var <- filter(
  exp1.Qangry.posterior.test,
  condition == c(
    "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
  )
) %>%
  droplevels()

par(mfrow = c(2, 2))

for(i in 1:length(levels(exp1.Qangry.posterior.test.var$condition))) {
  
  temp.data <- filter(exp1.Qangry.posterior.test.var, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

Take-home message:

• all faces are rated as less angry than angry faces, in both within-emotion variabilities
• angry and neutral faces are rated as angrier in high compared to low within-emotion variability
• happy and disgusted faces are rated as similarly angry in high and low within-emotion variability

---

---

How HAPPY are these faces?

# subset data
exp1.Qhappy.ratings <- ratings %>%
  filter(exp == "exp1", question == "happy") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
  ))

exp1.Qhappy.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp1.Qhappy.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
  adapt_delta = .99,
  max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)

saveRDS(exp1.Qhappy.model, file = "exp1_Qhappy_model.rds", compress = "gzip")

# subset data
exp1.Qhappy.ratings <- ratings %>%
  filter(exp == "exp1", question == "happy") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
  ))

# load model
exp1.Qhappy.model <- readRDS("exp1_Qhappy_model.rds")

# summary posterior distributions & diagnostics
summary.exp1.Qhappy.model <-
  describe_posterior(exp1.Qhappy.model,
    centrality = "mean",
    ci = 0.95,
    ci_method = "hdi",
    test = NULL,
    diagnostic = "all"
  ) %>%
  as_tibble() %>%
  separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
  dplyr::select(-coef) %>%
  mutate(Parameter = recode(factor(Parameter),
    "Intercept.1." = "Intercept[1]",
    "Intercept.2." = "Intercept[2]",
    "Intercept.3." = "Intercept[3]",
    "Intercept.4." = "Intercept[4]",
    "Intercept.5." = "Intercept[5]",
    "Intercept.6." = "Intercept[6]",
    "Intercept.7." = "Intercept[7]",
    "Intercept.8." = "Intercept[8]",
    "emotionanger.variabilityhigh" = "anger, high variability",
    "emotionneutral.variabilityhigh" = "neutral, high variability",
    "emotionhappiness.variabilityhigh" = "happiness, high variability",
    "emotiondisgust.variabilityhigh" = "disgust, high variability",
    "emotionanger.variabilitylow" = "anger, low variability",
    "emotionneutral.variabilitylow" = "neutral, low variability",
    "emotionhappiness.variabilitylow" = "happiness, low variability",
    "emotiondisgust.variabilitylow" = "disgust, low variability"
  )) %>%
  dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)

kable(summary.exp1.Qhappy.model, digits = 2)

Parameter	Mean	HDI95_Low	HDI95_High	Estimated_Sample_Size	MonteCarlo_StdErr
Intercept[1]	-1.41	-3.60	0.78	6491	0.01
Intercept[2]	-0.41	-2.58	1.79	6516	0.01
Intercept[3]	0.48	-1.80	2.68	6386	0.01
Intercept[4]	0.84	-1.35	3.14	6338	0.01
Intercept[5]	2.01	-0.28	4.33	6408	0.01
Intercept[6]	2.46	0.17	4.83	6484	0.01
Intercept[7]	4.65	2.17	7.14	6622	0.02
Intercept[8]	6.58	3.74	9.46	6104	0.02
happiness, high variability	5.75	3.11	8.46	6326	0.02
neutral, high variability	-1.63	-3.92	0.55	6560	0.01
anger, high variability	-2.41	-4.63	-0.18	6557	0.01
disgust, high variability	-2.91	-5.27	-0.54	6647	0.01
happiness, low variability	5.79	3.15	8.49	6596	0.02
neutral, low variability	-0.69	-2.89	1.53	6549	0.01
anger, low variability	-1.81	-4.08	0.39	6539	0.01
disgust, low variability	-2.25	-4.62	0.08	6649	0.01

# MCMC chains
exp1.Qhappy.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select("Chain" = ".chain",
                "Intercept[1]" = "b_Intercept[1]",
                "Intercept[2]" = "b_Intercept[2]",
                "Intercept[3]" = "b_Intercept[3]",
                "Intercept[4]" = "b_Intercept[4]",
                "Intercept[5]" = "b_Intercept[5]",
                "Intercept[6]" = "b_Intercept[6]",
                "Intercept[7]" = "b_Intercept[7]",
                "Intercept[8]" = "b_Intercept[8]",
                "anger, high variability" = "b_emotionanger:variabilityhigh",
                "neutral, high variability" = "b_emotionneutral:variabilityhigh",
                "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
                "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
                "anger, low variability" = "b_emotionanger:variabilitylow",
                "neutral, low variability" = "b_emotionneutral:variabilitylow",
                "happiness, low variability" = "b_emotionhappiness:variabilitylow",
                "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  as.data.frame() %>%
  mcmc_trace(facet_args = list(ncol = 2)) +
  geom_vline(xintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5
  ) +
  ggtitle("Chain convergence") +
  theme_EmoSSR +
  theme(legend.position = "none")

### observed vs. predicted data ###
exp1.Qhappy.ratings %>%
  dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
  mutate(predicted = posterior_predict(exp1.Qhappy.model)[1, ]) %>% # take only first row of sampled ratings (as example)
  dplyr::rename(observed = ratings) %>%
  gather(key = data_type, value = ratings, observed:predicted) %>%
  ggplot(aes(x = condition, 
             y = ratings,
             color = condition,
             fill = condition)) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(limits = c(1, 9),
                     breaks = seq(1, 9, 1)) +
  geom_hline(yintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5) +
  facet_wrap( ~ data_type, scales = "free") +
  xlab("") +
  ggtitle("experiment 1, \"happiness\" question ") +
  theme_EmoSSR

exp1.Qhappy.posterior.test <- exp1.Qhappy.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select(
    "anger, high variability" = "b_emotionanger:variabilityhigh",
    "neutral, high variability" = "b_emotionneutral:variabilityhigh",
    "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
    "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
    "anger, low variability" = "b_emotionanger:variabilitylow",
    "neutral, low variability" = "b_emotionneutral:variabilitylow",
    "happiness, low variability" = "b_emotionhappiness:variabilitylow",
    "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  mutate(
    ### emotion (happiness minus others)
    ## high variability
    # happiness minus neutral
    diff.high.happinessVSneutral = `happiness, high variability` - `neutral, high variability`,
    # happiness minus anger
    diff.high.happinessVSanger = `happiness, high variability` - `anger, high variability`,
    # happiness minus disgust
    diff.high.happinessVSdisgust = `happiness, high variability` - `disgust, high variability`,
    ## low variability
    # happiness minus neutral
    diff.low.happinessVSneutral = `happiness, low variability` - `neutral, low variability`,
    # happiness minus anger
    diff.low.happinessVSanger = `happiness, low variability` - `anger, low variability`,
    # happiness minus disgust
    diff.low.happinessVSdisgust = `happiness, low variability` - `disgust, low variability`,
    ### variability (high minus low)
    # anger
    diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
    # neutral
    diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
    # happiness
    diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
    # disgust
    diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
  ) %>%
  gather(key = "condition", value = "ratings") %>%
  filter(condition == c(
    "diff.high.happinessVSneutral",
    "diff.high.happinessVSanger",
    "diff.high.happinessVSdisgust",
    "diff.low.happinessVSneutral",
    "diff.low.happinessVSanger",
    "diff.low.happinessVSdisgust",
    "diff.anger.highVSlow",
    "diff.neutral.highVSlow",
    "diff.happiness.highVSlow",
    "diff.disgust.highVSlow"
  )) %>%
  mutate(condition = factor(condition,
    levels = c(
      "diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
      "diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust",
      "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
    )
  ))

# summary
summary.exp1.Qhappy.posterior.test <- exp1.Qhappy.posterior.test %>%
  group_by(condition) %>%
  dplyr::summarize(
    mean = mean(ratings),
    hdi.95.low = hdi(ratings)[1],
    hdi.95.high = hdi(ratings)[2],
    evidence.ratio = ifelse(mean < 0, 
                            length(which(ratings < 0)) / length(which(ratings > 0)),
                            length(which(ratings > 0)) / length(which(ratings < 0))
    )
  ) %>%
  ungroup()

kable(summary.exp1.Qhappy.posterior.test, digits = 2)

condition	mean	hdi.95.low	hdi.95.high	evidence.ratio
diff.high.happinessVSneutral	7.41	5.42	9.61	Inf
diff.high.happinessVSanger	8.15	5.89	10.20	Inf
diff.high.happinessVSdisgust	8.62	6.38	11.02	Inf
diff.low.happinessVSneutral	6.50	4.47	8.68	Inf
diff.low.happinessVSanger	7.60	5.67	9.93	Inf
diff.low.happinessVSdisgust	8.03	5.82	10.44	Inf
diff.anger.highVSlow	-0.61	-1.33	0.23	14.84
diff.neutral.highVSlow	-0.93	-1.79	-0.12	176.78
diff.happiness.highVSlow	-0.03	-1.28	0.97	1.07
diff.disgust.highVSlow	-0.68	-2.08	0.79	6.17

### Happiness vs. others ###

exp1.Qhappy.posterior.test.emoref <- filter(
  exp1.Qhappy.posterior.test,
  condition == c(
    "diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
    "diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust"
  )
) %>%
  droplevels()

par(mfrow = c(3, 2))

for(i in 1:length(levels(exp1.Qhappy.posterior.test.emoref$condition))) {
  
  temp.data <- filter(exp1.Qhappy.posterior.test.emoref, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

### High vs. low variability ###

exp1.Qhappy.posterior.test.var <- filter(
  exp1.Qhappy.posterior.test,
  condition == c(
    "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
  )
) %>%
  droplevels()

par(mfrow = c(2, 2))

for(i in 1:length(levels(exp1.Qhappy.posterior.test.var$condition))) {
  
  temp.data <- filter(exp1.Qhappy.posterior.test.var, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

Take-home message:

• all faces are rated as less happy than happy faces, in both within-emotion variabilities
• angry and neutral faces are rated as happier in low compared to high within-emotion variability
• happy and disgusted faces are rated as similarly happy in high and low within-emotion variability

---

---

How DISGUSTED are these faces?

# subset data
exp1.Qdisgust.ratings <- ratings %>%
  filter(exp == "exp1", question == "disgusted") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
  ))

exp1.Qdisgust.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp1.Qdisgust.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
  adapt_delta = .99,
  max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)

saveRDS(exp1.Qdisgust.model, file = "exp1_Qdisgust_model.rds", compress = "gzip")

# subset data
exp1.Qdisgust.ratings <- ratings %>%
  filter(exp == "exp1", question == "disgusted") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
  ))

# load model
exp1.Qdisgust.model <- readRDS("exp1_Qdisgust_model.rds")

# summary posterior distributions & diagnostics
summary.exp1.Qdisgust.model <-
  describe_posterior(exp1.Qdisgust.model,
    centrality = "mean",
    ci = 0.95,
    ci_method = "hdi",
    test = NULL,
    diagnostic = "all"
  ) %>%
  as_tibble() %>%
  separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
  dplyr::select(-coef) %>%
  mutate(Parameter = recode(factor(Parameter),
    "Intercept.1." = "Intercept[1]",
    "Intercept.2." = "Intercept[2]",
    "Intercept.3." = "Intercept[3]",
    "Intercept.4." = "Intercept[4]",
    "Intercept.5." = "Intercept[5]",
    "Intercept.6." = "Intercept[6]",
    "Intercept.7." = "Intercept[7]",
    "Intercept.8." = "Intercept[8]",
    "emotionanger.variabilityhigh" = "anger, high variability",
    "emotionneutral.variabilityhigh" = "neutral, high variability",
    "emotionhappiness.variabilityhigh" = "happiness, high variability",
    "emotiondisgust.variabilityhigh" = "disgust, high variability",
    "emotionanger.variabilitylow" = "anger, low variability",
    "emotionneutral.variabilitylow" = "neutral, low variability",
    "emotionhappiness.variabilitylow" = "happiness, low variability",
    "emotiondisgust.variabilitylow" = "disgust, low variability"
  )) %>%
  dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)

kable(summary.exp1.Qdisgust.model, digits = 2)

Parameter	Mean	HDI95_Low	HDI95_High	Estimated_Sample_Size	MonteCarlo_StdErr
Intercept[1]	-1.11	-3.26	1.02	7706	0.01
Intercept[2]	-0.41	-2.57	1.68	7889	0.01
Intercept[3]	0.14	-1.98	2.27	7925	0.01
Intercept[4]	0.70	-1.43	2.82	7905	0.01
Intercept[5]	1.36	-0.82	3.46	7833	0.01
Intercept[6]	1.90	-0.22	4.11	7703	0.01
Intercept[7]	2.85	0.56	5.00	7245	0.01
Intercept[8]	4.42	2.08	6.81	6244	0.02
disgust, high variability	3.83	1.38	6.18	7161	0.01
neutral, high variability	-0.47	-2.60	1.68	8064	0.01
anger, high variability	-0.17	-2.33	1.89	8057	0.01
happiness, high variability	-2.48	-5.11	0.00	8998	0.01
disgust, low variability	3.16	0.82	5.51	7563	0.01
neutral, low variability	-0.94	-3.11	1.22	8027	0.01
anger, low variability	0.24	-1.89	2.38	8179	0.01
happiness, low variability	-3.22	-5.83	-0.59	9266	0.01

# MCMC chains
exp1.Qdisgust.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select("Chain" = ".chain",
                "Intercept[1]" = "b_Intercept[1]",
                "Intercept[2]" = "b_Intercept[2]",
                "Intercept[3]" = "b_Intercept[3]",
                "Intercept[4]" = "b_Intercept[4]",
                "Intercept[5]" = "b_Intercept[5]",
                "Intercept[6]" = "b_Intercept[6]",
                "Intercept[7]" = "b_Intercept[7]",
                "Intercept[8]" = "b_Intercept[8]",
                "anger, high variability" = "b_emotionanger:variabilityhigh",
                "neutral, high variability" = "b_emotionneutral:variabilityhigh",
                "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
                "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
                "anger, low variability" = "b_emotionanger:variabilitylow",
                "neutral, low variability" = "b_emotionneutral:variabilitylow",
                "happiness, low variability" = "b_emotionhappiness:variabilitylow",
                "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  as.data.frame() %>%
  mcmc_trace(facet_args = list(ncol = 2)) +
  geom_vline(xintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5
  ) +
  ggtitle("Chain convergence") +
  theme_EmoSSR +
  theme(legend.position = "none")

### observed vs. predicted data ###
exp1.Qdisgust.ratings %>%
  dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
  mutate(predicted = posterior_predict(exp1.Qdisgust.model)[1, ]) %>% # take only first row of sampled ratings (as example)
  dplyr::rename(observed = ratings) %>%
  gather(key = data_type, value = ratings, observed:predicted) %>%
  ggplot(aes(x = condition, 
             y = ratings,
             color = condition,
             fill = condition)) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(limits = c(1, 9),
                     breaks = seq(1, 9, 1)) +
  geom_hline(yintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5) +
  facet_wrap( ~ data_type, scales = "free") +
  xlab("") +
  ggtitle("experiment 1, \"disgusted\" question ") +
  theme_EmoSSR

exp1.Qdisgust.posterior.test <- exp1.Qdisgust.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select(
    "anger, high variability" = "b_emotionanger:variabilityhigh",
    "neutral, high variability" = "b_emotionneutral:variabilityhigh",
    "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
    "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
    "anger, low variability" = "b_emotionanger:variabilitylow",
    "neutral, low variability" = "b_emotionneutral:variabilitylow",
    "happiness, low variability" = "b_emotionhappiness:variabilitylow",
    "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  mutate(
    ### emotion (disgust minus others)
    ## high variability
    # disgust minus neutral
    diff.high.disgustVSneutral = `disgust, high variability` - `neutral, high variability`,
    # disgust minus anger
    diff.high.disgustVSanger = `disgust, high variability` - `anger, high variability`,
    # disgust minus happiness
    diff.high.disgustVShappiness = `disgust, high variability` - `happiness, high variability`,
    ## low variability
    # disgust minus neutral
    diff.low.disgustVSneutral = `disgust, low variability` - `neutral, low variability`,
    # disgust minus anger
    diff.low.disgustVSanger = `disgust, low variability` - `anger, low variability`,
    # disgust minus happiness
    diff.low.disgustVShappiness = `disgust, low variability` - `happiness, low variability`,
    ### variability (high minus low)
    # anger
    diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
    # neutral
    diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
    # happiness
    diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
    # disgust
    diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
  ) %>%
  gather(key = "condition", value = "ratings") %>%
  filter(condition == c(
    "diff.high.disgustVSneutral",
    "diff.high.disgustVSanger",
    "diff.high.disgustVShappiness",
    "diff.low.disgustVSneutral",
    "diff.low.disgustVSanger",
    "diff.low.disgustVShappiness",
    "diff.anger.highVSlow",
    "diff.neutral.highVSlow",
    "diff.happiness.highVSlow",
    "diff.disgust.highVSlow"
  )) %>%
  mutate(condition = factor(condition,
    levels = c(
      "diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
      "diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness",
      "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
    )
  ))

# summary
summary.exp1.Qdisgust.posterior.test <- exp1.Qdisgust.posterior.test %>%
  group_by(condition) %>%
  dplyr::summarize(
    mean = mean(ratings),
    hdi.95.low = hdi(ratings)[1],
    hdi.95.high = hdi(ratings)[2],
    evidence.ratio = ifelse(mean < 0, 
                            length(which(ratings < 0)) / length(which(ratings > 0)),
                            length(which(ratings > 0)) / length(which(ratings < 0))
    )
  ) %>%
  ungroup()

kable(summary.exp1.Qdisgust.posterior.test, digits = 2)

condition	mean	hdi.95.low	hdi.95.high	evidence.ratio
diff.high.disgustVSneutral	4.31	2.94	5.98	Inf
diff.high.disgustVSanger	3.99	2.73	5.45	Inf
diff.high.disgustVShappiness	6.31	3.97	8.79	Inf
diff.low.disgustVSneutral	4.11	2.71	5.57	Inf
diff.low.disgustVSanger	2.93	1.66	4.17	Inf
diff.low.disgustVShappiness	6.37	4.15	8.85	Inf
diff.anger.highVSlow	-0.41	-1.07	0.25	7.60
diff.neutral.highVSlow	0.48	-0.34	1.27	7.94
diff.happiness.highVSlow	0.73	-1.46	2.74	3.32
diff.disgust.highVSlow	0.68	-0.37	1.96	8.14

### Disgust vs. others ###

exp1.Qdisgust.posterior.test.emoref <- filter(
  exp1.Qdisgust.posterior.test,
  condition == c(
    "diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
    "diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness"
  )
) %>%
  droplevels()

par(mfrow = c(3, 2))

for(i in 1:length(levels(exp1.Qdisgust.posterior.test.emoref$condition))) {
  
  temp.data <- filter(exp1.Qdisgust.posterior.test.emoref, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

### High vs. low variability ###

exp1.Qdisgust.posterior.test.var <- filter(
  exp1.Qdisgust.posterior.test,
  condition == c(
    "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
  )
) %>%
  droplevels()

par(mfrow = c(2, 2))

for(i in 1:length(levels(exp1.Qdisgust.posterior.test.var$condition))) {
  
  temp.data <- filter(exp1.Qdisgust.posterior.test.var, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

Take-home message:

• all faces are rated as less disgusted than disgusted faces, in both within-emotion variabilities
• all faces are rated as similarly disgusted in high and low within-emotion variability

---

---

Experiment 2

• 2 face identities
• 3 emotions (anger, happiness, disgust)

How ANGRY are these faces?

# subset data
exp2.Qangry.ratings <- ratings %>%
  filter(exp == "exp2", question == "angry") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
  ))

exp2.Qangry.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp2.Qangry.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
  adapt_delta = .99,
  max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)

saveRDS(exp2.Qangry.model, file = "exp2_Qangry_model.rds", compress = "gzip")

# subset data
exp2.Qangry.ratings <- ratings %>%
  filter(exp == "exp2", question == "angry") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
  ))

# load model
exp2.Qangry.model <- readRDS("exp2_Qangry_model.rds")

# summary posterior distributions & diagnostics
summary.exp2.Qangry.model <-
  describe_posterior(exp2.Qangry.model,
    centrality = "mean",
    ci = 0.95,
    ci_method = "hdi",
    test = NULL,
    diagnostic = "all"
  ) %>%
  as_tibble() %>%
  separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
  dplyr::select(-coef) %>%
  mutate(Parameter = recode(factor(Parameter),
    "Intercept.1." = "Intercept[1]",
    "Intercept.2." = "Intercept[2]",
    "Intercept.3." = "Intercept[3]",
    "Intercept.4." = "Intercept[4]",
    "Intercept.5." = "Intercept[5]",
    "Intercept.6." = "Intercept[6]",
    "Intercept.7." = "Intercept[7]",
    "Intercept.8." = "Intercept[8]",
    "emotionanger.variabilityhigh" = "anger, high variability",
    "emotionneutral.variabilityhigh" = "neutral, high variability",
    "emotionhappiness.variabilityhigh" = "happiness, high variability",
    "emotiondisgust.variabilityhigh" = "disgust, high variability",
    "emotionanger.variabilitylow" = "anger, low variability",
    "emotionneutral.variabilitylow" = "neutral, low variability",
    "emotionhappiness.variabilitylow" = "happiness, low variability",
    "emotiondisgust.variabilitylow" = "disgust, low variability"
  )) %>%
  dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)

kable(summary.exp2.Qangry.model, digits = 2)

Parameter	Mean	HDI95_Low	HDI95_High	Estimated_Sample_Size	MonteCarlo_StdErr
Intercept[1]	-0.85	-2.94	1.34	6243	0.01
Intercept[2]	-0.12	-2.26	2.02	6362	0.01
Intercept[3]	0.44	-1.66	2.62	6361	0.01
Intercept[4]	0.90	-1.20	3.08	6286	0.01
Intercept[5]	1.77	-0.33	4.01	6081	0.01
Intercept[6]	2.45	0.25	4.62	5836	0.01
Intercept[7]	3.32	1.07	5.53	5711	0.02
Intercept[8]	4.28	2.09	6.70	5549	0.02
anger, high variability	4.56	2.25	6.84	5539	0.02
neutral, high variability	0.69	-1.50	2.83	6383	0.01
happiness, high variability	-3.70	-6.66	-0.79	9174	0.02
disgust, high variability	0.42	-1.90	2.51	6832	0.01
anger, low variability	2.39	0.21	4.56	5951	0.01
neutral, low variability	-0.36	-2.40	1.89	6419	0.01
happiness, low variability	-4.59	-7.96	-1.48	8720	0.02
disgust, low variability	0.43	-1.76	2.53	6519	0.01

# MCMC chains
exp2.Qangry.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select("Chain" = ".chain",
                "Intercept[1]" = "b_Intercept[1]",
                "Intercept[2]" = "b_Intercept[2]",
                "Intercept[3]" = "b_Intercept[3]",
                "Intercept[4]" = "b_Intercept[4]",
                "Intercept[5]" = "b_Intercept[5]",
                "Intercept[6]" = "b_Intercept[6]",
                "Intercept[7]" = "b_Intercept[7]",
                "Intercept[8]" = "b_Intercept[8]",
                "anger, high variability" = "b_emotionanger:variabilityhigh",
                "neutral, high variability" = "b_emotionneutral:variabilityhigh",
                "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
                "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
                "anger, low variability" = "b_emotionanger:variabilitylow",
                "neutral, low variability" = "b_emotionneutral:variabilitylow",
                "happiness, low variability" = "b_emotionhappiness:variabilitylow",
                "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  as.data.frame() %>%
  mcmc_trace(facet_args = list(ncol = 2)) +
  geom_vline(xintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5
  ) +
  ggtitle("Chain convergence") +
  theme_EmoSSR +
  theme(legend.position = "none")

### observed vs. predicted data ###
exp2.Qangry.ratings %>%
  dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
  mutate(predicted = posterior_predict(exp2.Qangry.model)[1, ]) %>% # take only first row of sampled ratings (as example)
  dplyr::rename(observed = ratings) %>%
  gather(key = data_type, value = ratings, observed:predicted) %>%
  ggplot(aes(x = condition, 
             y = ratings,
             color = condition,
             fill = condition)) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(limits = c(1, 9),
                     breaks = seq(1, 9, 1)) +
  geom_hline(yintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5) +
  facet_wrap( ~ data_type, scales = "free") +
  xlab("") +
  ggtitle("pilot 2, \"anger\" question ") +
  theme_EmoSSR

exp2.Qangry.posterior.test <- exp2.Qangry.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select(
    "anger, high variability" = "b_emotionanger:variabilityhigh",
    "neutral, high variability" = "b_emotionneutral:variabilityhigh",
    "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
    "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
    "anger, low variability" = "b_emotionanger:variabilitylow",
    "neutral, low variability" = "b_emotionneutral:variabilitylow",
    "happiness, low variability" = "b_emotionhappiness:variabilitylow",
    "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  mutate(
    ### emotion (anger minus others)
    ## high variability
    # anger minus neutral
    diff.high.angerVSneutral = `anger, high variability` - `neutral, high variability`,
    # anger minus happiness
    diff.high.angerVShappiness = `anger, high variability` - `happiness, high variability`,
    # anger minus disgust
    diff.high.angerVSdisgust = `anger, high variability` - `disgust, high variability`,
    ## low variability
    # anger minus neutral
    diff.low.angerVSneutral = `anger, low variability` - `neutral, low variability`,
    # anger minus happiness
    diff.low.angerVShappiness = `anger, low variability` - `happiness, low variability`,
    # anger minus disgust
    diff.low.angerVSdisgust = `anger, low variability` - `disgust, low variability`,
    ### variability (high minus low)
    # anger
    diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
    # neutral
    diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
    # happiness
    diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
    # disgust
    diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
  ) %>%
  gather(key = "condition", value = "ratings") %>%
  filter(condition == c(
    "diff.high.angerVSneutral",
    "diff.high.angerVShappiness",
    "diff.high.angerVSdisgust",
    "diff.low.angerVSneutral",
    "diff.low.angerVShappiness",
    "diff.low.angerVSdisgust",
    "diff.anger.highVSlow",
    "diff.neutral.highVSlow",
    "diff.happiness.highVSlow",
    "diff.disgust.highVSlow"
  )) %>%
  mutate(condition = factor(condition,
    levels = c(
      "diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
      "diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust",
      "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
    )
  ))

# summary
summary.exp2.Qangry.posterior.test <- exp2.Qangry.posterior.test %>%
  group_by(condition) %>%
  dplyr::summarize(
    mean = mean(ratings),
    hdi.95.low = hdi(ratings)[1],
    hdi.95.high = hdi(ratings)[2],
    evidence.ratio = ifelse(mean < 0, 
                            length(which(ratings < 0)) / length(which(ratings > 0)),
                            length(which(ratings > 0)) / length(which(ratings < 0))
    )
  ) %>%
  ungroup()

kable(summary.exp2.Qangry.posterior.test, digits = 2)

condition	mean	hdi.95.low	hdi.95.high	evidence.ratio
diff.high.angerVSneutral	3.87	2.88	4.96	Inf
diff.high.angerVShappiness	8.22	6.03	11.20	Inf
diff.high.angerVSdisgust	4.15	2.86	5.44	Inf
diff.low.angerVSneutral	2.77	1.88	3.63	Inf
diff.low.angerVShappiness	6.94	4.19	9.75	Inf
diff.low.angerVSdisgust	1.96	1.21	2.85	Inf
diff.anger.highVSlow	2.17	1.47	2.99	Inf
diff.neutral.highVSlow	1.04	0.25	1.76	319.00
diff.happiness.highVSlow	0.87	-2.62	4.57	2.08
diff.disgust.highVSlow	-0.01	-0.96	0.95	1.04

### Anger vs. others ###

exp2.Qangry.posterior.test.emoref <- filter(
  exp2.Qangry.posterior.test,
  condition == c(
    "diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
    "diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust"
  )
) %>%
  droplevels()

par(mfrow = c(3, 2))

for(i in 1:length(levels(exp2.Qangry.posterior.test.emoref$condition))) {
  
  temp.data <- filter(exp2.Qangry.posterior.test.emoref, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

### High vs. low variability ###

exp2.Qangry.posterior.test.var <- filter(
  exp2.Qangry.posterior.test,
  condition == c(
    "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
  )
) %>%
  droplevels()

par(mfrow = c(2, 2))

for(i in 1:length(levels(exp2.Qangry.posterior.test.var$condition))) {
  
  temp.data <- filter(exp2.Qangry.posterior.test.var, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

Take-home message:

• all faces are rated as less angry than angry faces, in both within-emotion variabilities
• angry and neutral faces are rated as angrier in high compared to low within-emotion variability
• happy and disgusted faces are rated as similarly angry in high and low within-emotion variability

---

---

How HAPPY are these faces?

# subset data
exp2.Qhappy.ratings <- ratings %>%
  filter(exp == "exp2", question == "happy") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
  ))

exp2.Qhappy.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp2.Qhappy.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
  adapt_delta = .99,
  max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)

saveRDS(exp2.Qhappy.model, file = "exp2_Qhappy_model.rds", compress = "gzip")

# subset data
exp2.Qhappy.ratings <- ratings %>%
  filter(exp == "exp2", question == "happy") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
  ))

# load model
exp2.Qhappy.model <- readRDS("exp2_Qhappy_model.rds")

# summary posterior distributions & diagnostics
summary.exp2.Qhappy.model <-
  describe_posterior(exp2.Qhappy.model,
    centrality = "mean",
    ci = 0.95,
    ci_method = "hdi",
    test = NULL,
    diagnostic = "all"
  ) %>%
  as_tibble() %>%
  separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
  dplyr::select(-coef) %>%
  mutate(Parameter = recode(factor(Parameter),
    "Intercept.1." = "Intercept[1]",
    "Intercept.2." = "Intercept[2]",
    "Intercept.3." = "Intercept[3]",
    "Intercept.4." = "Intercept[4]",
    "Intercept.5." = "Intercept[5]",
    "Intercept.6." = "Intercept[6]",
    "Intercept.7." = "Intercept[7]",
    "Intercept.8." = "Intercept[8]",
    "emotionanger.variabilityhigh" = "anger, high variability",
    "emotionneutral.variabilityhigh" = "neutral, high variability",
    "emotionhappiness.variabilityhigh" = "happiness, high variability",
    "emotiondisgust.variabilityhigh" = "disgust, high variability",
    "emotionanger.variabilitylow" = "anger, low variability",
    "emotionneutral.variabilitylow" = "neutral, low variability",
    "emotionhappiness.variabilitylow" = "happiness, low variability",
    "emotiondisgust.variabilitylow" = "disgust, low variability"
  )) %>%
  dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)

kable(summary.exp2.Qhappy.model, digits = 2)

Parameter	Mean	HDI95_Low	HDI95_High	Estimated_Sample_Size	MonteCarlo_StdErr
Intercept[1]	-1.51	-3.73	0.65	5492	0.02
Intercept[2]	-0.57	-2.74	1.62	5637	0.01
Intercept[3]	0.26	-1.91	2.43	5720	0.01
Intercept[4]	0.65	-1.51	2.86	5879	0.01
Intercept[5]	2.12	-0.14	4.30	5996	0.01
Intercept[6]	2.74	0.52	5.04	6054	0.01
Intercept[7]	3.95	1.58	6.32	5886	0.02
Intercept[8]	5.85	3.26	8.57	5132	0.02
happiness, high variability	5.55	2.98	8.19	5577	0.02
neutral, high variability	-0.67	-2.87	1.65	5954	0.02
anger, high variability	-2.71	-5.01	-0.57	5506	0.02
disgust, high variability	-3.02	-5.26	-0.66	5617	0.02
happiness, low variability	5.08	2.46	7.62	5763	0.02
neutral, low variability	0.19	-1.96	2.36	5764	0.01
anger, low variability	-1.78	-4.00	0.41	5573	0.02
disgust, low variability	-2.75	-5.01	-0.47	5530	0.02

# MCMC chains
exp2.Qhappy.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select("Chain" = ".chain",
                "Intercept[1]" = "b_Intercept[1]",
                "Intercept[2]" = "b_Intercept[2]",
                "Intercept[3]" = "b_Intercept[3]",
                "Intercept[4]" = "b_Intercept[4]",
                "Intercept[5]" = "b_Intercept[5]",
                "Intercept[6]" = "b_Intercept[6]",
                "Intercept[7]" = "b_Intercept[7]",
                "Intercept[8]" = "b_Intercept[8]",
                "anger, high variability" = "b_emotionanger:variabilityhigh",
                "neutral, high variability" = "b_emotionneutral:variabilityhigh",
                "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
                "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
                "anger, low variability" = "b_emotionanger:variabilitylow",
                "neutral, low variability" = "b_emotionneutral:variabilitylow",
                "happiness, low variability" = "b_emotionhappiness:variabilitylow",
                "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  as.data.frame() %>%
  mcmc_trace(facet_args = list(ncol = 2)) +
  geom_vline(xintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5
  ) +
  ggtitle("Chain convergence") +
  theme_EmoSSR +
  theme(legend.position = "none")

### observed vs. predicted data ###
exp2.Qhappy.ratings %>%
  dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
  mutate(predicted = posterior_predict(exp2.Qhappy.model)[1, ]) %>% # take only first row of sampled ratings (as example)
  dplyr::rename(observed = ratings) %>%
  gather(key = data_type, value = ratings, observed:predicted) %>%
  ggplot(aes(x = condition, 
             y = ratings,
             color = condition,
             fill = condition)) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(limits = c(1, 9),
                     breaks = seq(1, 9, 1)) +
  geom_hline(yintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5) +
  facet_wrap( ~ data_type, scales = "free") +
  xlab("") +
  ggtitle("experiment 2, \"happiness\" question ") +
  theme_EmoSSR

exp2.Qhappy.posterior.test <- exp2.Qhappy.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select(
    "anger, high variability" = "b_emotionanger:variabilityhigh",
    "neutral, high variability" = "b_emotionneutral:variabilityhigh",
    "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
    "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
    "anger, low variability" = "b_emotionanger:variabilitylow",
    "neutral, low variability" = "b_emotionneutral:variabilitylow",
    "happiness, low variability" = "b_emotionhappiness:variabilitylow",
    "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  mutate(
    ### emotion (happiness minus others)
    ## high variability
    # happiness minus neutral
    diff.high.happinessVSneutral = `happiness, high variability` - `neutral, high variability`,
    # happiness minus anger
    diff.high.happinessVSanger = `happiness, high variability` - `anger, high variability`,
    # happiness minus disgust
    diff.high.happinessVSdisgust = `happiness, high variability` - `disgust, high variability`,
    ## low variability
    # happiness minus neutral
    diff.low.happinessVSneutral = `happiness, low variability` - `neutral, low variability`,
    # happiness minus anger
    diff.low.happinessVSanger = `happiness, low variability` - `anger, low variability`,
    # happiness minus disgust
    diff.low.happinessVSdisgust = `happiness, low variability` - `disgust, low variability`,
    ### variability (high minus low)
    # anger
    diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
    # neutral
    diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
    # happiness
    diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
    # disgust
    diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
  ) %>%
  gather(key = "condition", value = "ratings") %>%
  filter(condition == c(
    "diff.high.happinessVSneutral",
    "diff.high.happinessVSanger",
    "diff.high.happinessVSdisgust",
    "diff.low.happinessVSneutral",
    "diff.low.happinessVSanger",
    "diff.low.happinessVSdisgust",
    "diff.anger.highVSlow",
    "diff.neutral.highVSlow",
    "diff.happiness.highVSlow",
    "diff.disgust.highVSlow"
  )) %>%
  mutate(condition = factor(condition,
    levels = c(
      "diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
      "diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust",
      "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
    )
  ))

# summary
summary.exp2.Qhappy.posterior.test <- exp2.Qhappy.posterior.test %>%
  group_by(condition) %>%
  dplyr::summarize(
    mean = mean(ratings),
    hdi.95.low = hdi(ratings)[1],
    hdi.95.high = hdi(ratings)[2],
    evidence.ratio = ifelse(mean < 0, 
                            length(which(ratings < 0)) / length(which(ratings > 0)),
                            length(which(ratings > 0)) / length(which(ratings < 0))
    )
  ) %>%
  ungroup()

kable(summary.exp2.Qhappy.posterior.test, digits = 2)

condition	mean	hdi.95.low	hdi.95.high	evidence.ratio
diff.high.happinessVSneutral	6.23	4.18	8.51	Inf
diff.high.happinessVSanger	8.25	6.11	10.62	Inf
diff.high.happinessVSdisgust	8.55	6.35	11.02	Inf
diff.low.happinessVSneutral	4.90	3.09	6.75	Inf
diff.low.happinessVSanger	6.88	4.81	8.79	Inf
diff.low.happinessVSdisgust	7.81	5.59	9.87	Inf
diff.anger.highVSlow	-0.93	-1.70	-0.07	71.73
diff.neutral.highVSlow	-0.85	-1.74	0.29	24.00
diff.happiness.highVSlow	0.46	-0.99	1.86	3.07
diff.disgust.highVSlow	-0.27	-1.55	0.77	2.26

### Happiness vs. others ###

exp2.Qhappy.posterior.test.emoref <- filter(
  exp2.Qhappy.posterior.test,
  condition == c(
    "diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
    "diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust"
  )
) %>%
  droplevels()

par(mfrow = c(3, 2))

for(i in 1:length(levels(exp2.Qhappy.posterior.test.emoref$condition))) {
  
  temp.data <- filter(exp2.Qhappy.posterior.test.emoref, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

### High vs. low variability ###

exp2.Qhappy.posterior.test.var <- filter(
  exp2.Qhappy.posterior.test,
  condition == c(
    "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
  )
) %>%
  droplevels()

par(mfrow = c(2, 2))

for(i in 1:length(levels(exp2.Qhappy.posterior.test.var$condition))) {
  
  temp.data <- filter(exp2.Qhappy.posterior.test.var, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

Take-home message:

• all faces are rated as less happy than happy faces, in both within-emotion variabilities
• angry and neutral faces are rated as happier in low compared to high within-emotion variability
• happy and disgusted faces are rated as similarly happy in high and low within-emotion variability

---

---

How DISGUSTED are these faces?

# subset data
exp2.Qdisgust.ratings <- ratings %>%
  filter(exp == "exp2", question == "disgusted") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
  ))

exp2.Qdisgust.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp2.Qdisgust.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
  adapt_delta = .99,
  max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)

saveRDS(exp2.Qdisgust.model, file = "exp2_Qdisgust_model.rds", compress = "gzip")

# subset data
exp2.Qdisgust.ratings <- ratings %>%
  filter(exp == "exp2", question == "disgusted") %>%
  droplevels() %>%
  mutate(emotion = factor(emotion,
    levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
  ))

# load model
exp2.Qdisgust.model <- readRDS("exp2_Qdisgust_model.rds")

# summary posterior distributions & diagnostics
summary.exp2.Qdisgust.model <-
  describe_posterior(exp2.Qdisgust.model,
    centrality = "mean",
    ci = 0.95,
    ci_method = "hdi",
    test = NULL,
    diagnostic = "all"
  ) %>%
  as_tibble() %>%
  separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
  dplyr::select(-coef) %>%
  mutate(Parameter = recode(factor(Parameter),
    "Intercept.1." = "Intercept[1]",
    "Intercept.2." = "Intercept[2]",
    "Intercept.3." = "Intercept[3]",
    "Intercept.4." = "Intercept[4]",
    "Intercept.5." = "Intercept[5]",
    "Intercept.6." = "Intercept[6]",
    "Intercept.7." = "Intercept[7]",
    "Intercept.8." = "Intercept[8]",
    "emotionanger.variabilityhigh" = "anger, high variability",
    "emotionneutral.variabilityhigh" = "neutral, high variability",
    "emotionhappiness.variabilityhigh" = "happiness, high variability",
    "emotiondisgust.variabilityhigh" = "disgust, high variability",
    "emotionanger.variabilitylow" = "anger, low variability",
    "emotionneutral.variabilitylow" = "neutral, low variability",
    "emotionhappiness.variabilitylow" = "happiness, low variability",
    "emotiondisgust.variabilitylow" = "disgust, low variability"
  )) %>%
  dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)

kable(summary.exp2.Qdisgust.model, digits = 2)

Parameter	Mean	HDI95_Low	HDI95_High	Estimated_Sample_Size	MonteCarlo_StdErr
Intercept[1]	-1.44	-3.55	0.72	7608	0.01
Intercept[2]	-0.39	-2.48	1.78	7759	0.01
Intercept[3]	0.15	-2.04	2.24	7783	0.01
Intercept[4]	0.46	-1.59	2.68	7819	0.01
Intercept[5]	1.18	-0.95	3.34	7805	0.01
Intercept[6]	1.71	-0.35	3.98	7765	0.01
Intercept[7]	2.23	0.00	4.36	7664	0.01
Intercept[8]	3.51	1.15	5.71	6736	0.01
disgust, high variability	3.99	1.54	6.51	7636	0.01
neutral, high variability	-0.05	-2.18	2.16	8025	0.01
anger, high variability	-0.48	-2.67	1.60	7827	0.01
happiness, high variability	-3.61	-6.31	-0.79	9212	0.01
disgust, low variability	3.07	0.73	5.31	7593	0.01
neutral, low variability	-1.11	-3.27	1.06	7765	0.01
anger, low variability	0.49	-1.66	2.60	7876	0.01
happiness, low variability	-2.29	-4.48	-0.01	7681	0.01

# MCMC chains
exp2.Qdisgust.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select("Chain" = ".chain",
                "Intercept[1]" = "b_Intercept[1]",
                "Intercept[2]" = "b_Intercept[2]",
                "Intercept[3]" = "b_Intercept[3]",
                "Intercept[4]" = "b_Intercept[4]",
                "Intercept[5]" = "b_Intercept[5]",
                "Intercept[6]" = "b_Intercept[6]",
                "Intercept[7]" = "b_Intercept[7]",
                "Intercept[8]" = "b_Intercept[8]",
                "anger, high variability" = "b_emotionanger:variabilityhigh",
                "neutral, high variability" = "b_emotionneutral:variabilityhigh",
                "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
                "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
                "anger, low variability" = "b_emotionanger:variabilitylow",
                "neutral, low variability" = "b_emotionneutral:variabilitylow",
                "happiness, low variability" = "b_emotionhappiness:variabilitylow",
                "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  as.data.frame() %>%
  mcmc_trace(facet_args = list(ncol = 2)) +
  geom_vline(xintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5
  ) +
  ggtitle("Chain convergence") +
  theme_EmoSSR +
  theme(legend.position = "none")

### observed vs. predicted data ###
exp2.Qdisgust.ratings %>%
  dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
  mutate(predicted = posterior_predict(exp2.Qdisgust.model)[1, ]) %>% # take only first row of sampled ratings (as example)
  dplyr::rename(observed = ratings) %>%
  gather(key = data_type, value = ratings, observed:predicted) %>%
  ggplot(aes(x = condition, 
             y = ratings,
             color = condition,
             fill = condition)) +
  geom_pirate(
    bars = FALSE,
    cis = TRUE,
    lines = TRUE, lines_params = list(color = "black"),
    points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
    violins = TRUE, violins_params = list(size = 1),
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(option = "cividis") +
  scale_color_viridis_d(option = "cividis") +
  scale_y_continuous(limits = c(1, 9),
                     breaks = seq(1, 9, 1)) +
  geom_hline(yintercept = seq(1, 9, 1),
             linetype = "dotted",
             colour = "#999999",
             size = .8,
             alpha = .5) +
  facet_wrap( ~ data_type, scales = "free") +
  xlab("") +
  ggtitle("experiment 2, \"disgusted\" question ") +
  theme_EmoSSR

exp2.Qdisgust.posterior.test <- exp2.Qdisgust.model %>%
  spread_draws(`b_.*`, regex = TRUE) %>%
  dplyr::select(
    "anger, high variability" = "b_emotionanger:variabilityhigh",
    "neutral, high variability" = "b_emotionneutral:variabilityhigh",
    "happiness, high variability" = "b_emotionhappiness:variabilityhigh",
    "disgust, high variability" = "b_emotiondisgust:variabilityhigh",
    "anger, low variability" = "b_emotionanger:variabilitylow",
    "neutral, low variability" = "b_emotionneutral:variabilitylow",
    "happiness, low variability" = "b_emotionhappiness:variabilitylow",
    "disgust, low variability" = "b_emotiondisgust:variabilitylow"
  ) %>%
  mutate(
    ### emotion (disgust minus others)
    ## high variability
    # disgust minus neutral
    diff.high.disgustVSneutral = `disgust, high variability` - `neutral, high variability`,
    # disgust minus anger
    diff.high.disgustVSanger = `disgust, high variability` - `anger, high variability`,
    # disgust minus happiness
    diff.high.disgustVShappiness = `disgust, high variability` - `happiness, high variability`,
    ## low variability
    # disgust minus neutral
    diff.low.disgustVSneutral = `disgust, low variability` - `neutral, low variability`,
    # disgust minus anger
    diff.low.disgustVSanger = `disgust, low variability` - `anger, low variability`,
    # disgust minus happiness
    diff.low.disgustVShappiness = `disgust, low variability` - `happiness, low variability`,
    ### variability (high minus low)
    # anger
    diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
    # neutral
    diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
    # happiness
    diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
    # disgust
    diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
  ) %>%
  gather(key = "condition", value = "ratings") %>%
  filter(condition == c(
    "diff.high.disgustVSneutral",
    "diff.high.disgustVSanger",
    "diff.high.disgustVShappiness",
    "diff.low.disgustVSneutral",
    "diff.low.disgustVSanger",
    "diff.low.disgustVShappiness",
    "diff.anger.highVSlow",
    "diff.neutral.highVSlow",
    "diff.happiness.highVSlow",
    "diff.disgust.highVSlow"
  )) %>%
  mutate(condition = factor(condition,
    levels = c(
      "diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
      "diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness",
      "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
    )
  ))

# summary
summary.exp2.Qdisgust.posterior.test <- exp2.Qdisgust.posterior.test %>%
  group_by(condition) %>%
  dplyr::summarize(
    mean = mean(ratings),
    hdi.95.low = hdi(ratings)[1],
    hdi.95.high = hdi(ratings)[2],
    evidence.ratio = ifelse(mean < 0, 
                            length(which(ratings < 0)) / length(which(ratings > 0)),
                            length(which(ratings > 0)) / length(which(ratings < 0))
    )
  ) %>%
  ungroup()

kable(summary.exp2.Qdisgust.posterior.test, digits = 2)

condition	mean	hdi.95.low	hdi.95.high	evidence.ratio
diff.high.disgustVSneutral	4.02	2.32	5.57	Inf
diff.high.disgustVSanger	4.49	3.06	6.19	Inf
diff.high.disgustVShappiness	7.60	5.05	10.17	Inf
diff.low.disgustVSneutral	4.20	2.93	5.57	Inf
diff.low.disgustVSanger	2.60	1.52	3.60	Inf
diff.low.disgustVShappiness	5.38	3.92	6.94	Inf
diff.anger.highVSlow	-0.97	-1.56	-0.41	799.00
diff.neutral.highVSlow	1.07	0.19	1.86	159.00
diff.happiness.highVSlow	-1.30	-3.25	0.66	10.85
diff.disgust.highVSlow	0.93	-0.30	2.33	16.20

### Disgust vs. others ###

exp2.Qdisgust.posterior.test.emoref <- filter(
  exp2.Qdisgust.posterior.test,
  condition == c(
    "diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
    "diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness"
  )
) %>%
  droplevels()

par(mfrow = c(3, 2))

for(i in 1:length(levels(exp2.Qdisgust.posterior.test.emoref$condition))) {
  
  temp.data <- filter(exp2.Qdisgust.posterior.test.emoref, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

### High vs. low variability ###

exp2.Qdisgust.posterior.test.var <- filter(
  exp2.Qdisgust.posterior.test,
  condition == c(
    "diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
  )
) %>%
  droplevels()

par(mfrow = c(2, 2))

for(i in 1:length(levels(exp2.Qdisgust.posterior.test.var$condition))) {
  
  temp.data <- filter(exp2.Qdisgust.posterior.test.var, condition == levels(condition)[i])
  
  plotPost(
    pull(temp.data, ratings), 
    credMass = 0.95,
    compVal = 0,
    ROPE = NULL,
    xlab = "", 
    col = "#b3cde0", 
    showCurve = FALSE,
    cex = 1,
    main = levels(temp.data$condition)[i]
  )
  
}

Take-home message:

• all faces are rated as less disgusted than disgusted faces, in both within-emotion variabilities
• all faces are rated as similarly disgusted in high and low within-emotion variability

---

---